Categories

## Using Sympy’s Common Subexpression Elimination to generate Code

A nifty feature of sympy is its code generator. This allows generating code in many languages from expressions derived in sympy. I’ve been using this for a while and came across the Common Subexpression Elimination (CSE) function (again), recently.

In earlier* versions of sympy I did not find it of great use. But now I was pleasantly surprised how much it matured. (* I cannot say for sure when I last tried using it; it must’ve been a long while ago.)

What is it and what is it good for?

Let’s say we want to generate code for a PI controller in C. The literature says its transfer function (in DIN notation) reads:

Then furthermore the representation for a discrete time system with time step becomes

In order to efficiently and flexibly calculate the response we want to implement it in a second order IIR (infinite impulse response) filter structure. This structure is basically a representation of as a rational function

so we need to take the expression apart such that we can determine the coefficients a_1, a_2, b_0, … (by convention, the problem is scaled such that a_0 = 1.

This is all very cumbersome to do manually (and I suck at doing this, especially as I insert quite embarrassing mistakes all the time). So, sympy to the rescue:

import sympy as sy

Kp, Ki, Ts = sy.symbols("K_p K_i T_s", real=True)
z = sy.symbols("z", complex=True)
G_z = Kp + Ki*Ts/2 + (Ki*Ts/2 - Kp)/z

display(G_z)

n, d = G_z.ratsimp().as_numer_denom()
n, d = n.as_poly(z), d.as_poly(z)

display(n, d)

We first define the symbols needed. Then set up G_z according to our equation above. After displaying it for good measure, we do a rational simplification ratsimp() on it and split it into numerator n and denominator d.

Now, we represent n and d as polynomials of z. And this is how we obtain our coefficients a1, a2, and so on. We continue:

scale = d.LC()

n, d = (n.as_expr()/scale).as_poly(z), (d.as_expr()/scale).as_poly(z)
N, D = n.degree(), d.degree()

We find the scaling factor as the leading coefficient (LC()) of the denominator d, scale both accordingly and determine the degree of each (this is mostly to determine for how many terms we have to iterate when printing etc.).

Bs = n.all_coeffs()
As = d.all_coeffs()

display(Bs)
display(As)

Now, basically Bs contains all coefficients [b0 b1] and likewise for As. Here’s a shortcut, see if you can spot it.

And now comes the cse magic:

commons, exprs = sy.cse(Bs)

This returns a tuple. The first element is a list of tuples of arbitrarily assigned place holder variables x0, x1, and so on. The second element are the original expressions, but the placeholders are substituted. In this example:

(
[
(x0, K_i*T_s/2)
],
[K_p + x0, -K_p + x0]
)

Which means, instead of calculating K_i*T_s/2 twice, we can do it once, assign it to x0 and use that on the expressions.

You can even concatenate As and Bs for good measure (not really necessary here, since the As only contain 1 and 0, but in other controller topologies this may change).

commons, exprs = sy.cse(As + Bs)

# Printing C code for instance
for com in commons:
print( sy.ccode(com[1], assign_to=repr(com[0])) )

for expr, name in zip(exprs, "a0 a1 b0 b1".split(" ")):
print( sy.ccode(expr, assign_to=f"coeffs.{name:s}") )

We obtain the common terms commons and the reduced expressions exprs, then print the common terms assigning them to their respective placeholder variables, and finally print the reduced expressions. This produces this five liner in C:

x0 = (1.0/2.0)*K_i*T_s;
coeffs.a0 = 1;
coeffs.a1 = 0;
coeffs.b0 = K_p + x0;
coeffs.b1 = -K_p + x0;

Now, how is this useful you ask? Well, for this toy example it certainly is a lot of work for five lines I could’ve just typed out.

But if you go into a PID controller with damped differential term, and you want your code to be efficient (i.e. have as few operations as possible), this is very handy. Just by changing the symbols and the expression for G_z we obtain this:

x0 = 2*T_v;
x1 = 1.0/(T_s + x0);
x2 = 4*T_v;
x3 = T_i*T_s;
x4 = 2*x3;
x5 = T_i*x2;
x6 = 1.0/(x4 + x5);
x7 = K_p*x4;
x8 = K_p*T_s*x0;
x9 = pow(T_s, 2);
x10 = 4*K_p*T_d*T_i + K_p*x5;
x11 = K_p*x9 + x10;
c.a0 = 1;
c.a1 = -x1*x2;
c.a2 = x1*(-T_s + 2*T_v);
c.b0 = x6*(x11 + x7 + x8);
c.b1 = (K_p*x9 - x10)/(T_i*x0 + x3);
c.b2 = x6*(x11 - x7 - x8);

I did this manually once, and only obtained five or six of the placeholders. That’s the power of cse().

Oh, you need the code in Python for testing? No problem, just use pycode instead (and work around that assign_to parameter).

So, there you have it.

Categories

… for all the obvious reasons.

Since then, I have noticed that I had addictive behavioural patterns. Given the relative smallness (is that a word) of the Fediverse, I hope I will do more sensible stuff with my time.

I had noticed this before, when I ditched Facebook. So, bottom line: “social media” isn’t good for me.

Categories

## Got a R&S FSH4 and want to extract screenshots?

Situation: I’ve taken about 50ish measurements on an FSH4 spectrum analyser saving a dataset of each using the save data function of the spectrum analyser.

Turns out: The R&S InstrumentView application does not install on Wine. So, how do I get the screenshots embedded in the saved datasets? Right: Search for the PNG magic number and a chunk end.

This is really sketchy but worked on all measurement files. If you’re in the same situation, just check out my repository.

https://gitlab.com/cweickhmann/extract-png-from-rs-set

Note: No affiliation with R&S whatsoever. Please don’t pester them if the script is not working.

Categories

There are various possible ways you can lose access to your Gitea instance. In particular, if you are getting user accounts via LDAP, you should always have a key to get back in: a “local” admin.

In case you haven’t got one (totally not like I did definitely only once 🤫), deleted it or lost its password, here’s the procedure in to variants: running on plain Linux, and running within Docker.

## Plain Linux

List the available (local) user accounts (you may need to run this with sudo):

$gitea admin user list ID Username Email IsActive IsAdmin 1 user1 me@home.earth true false You may not get a single user listed here. In order to add a local administrator, do this $ sudo gitea admin user create --username local_admin --email admins@email.earth --admin --random-password

New user 'local_admin' has been successfully created!

$sudo gitea admin user list ID Username Email IsActive IsAdmin 1 unpriv_local bowfingermail@gmx.net true true 3 local_admin admins@email.earth true true Especially on shared machines, I’d strongly recommend not setting a password via the command line. It’ll stick in bash_history and may be visible in the process list. Hence the --random-password option here. Use the generated password upon first login in the browser. ## Docker Container In order to get to the “plain Linux” field, you simply run a bash in the Gitea Docker container. What’s the container’s name? Get it via docker container ls -a. Let’s say it’s called gitea. To start bash in that container, do user1@linux:$ docker exec -it gitea /bin/bash

Now, you basically do the same thing as in the section above.

docker # gitea admin user list

1    user1        me@home.earth     true     false

New user 'local_admin' has been successfully created!

docker # gitea admin user list

1    unpriv_local bowfingermail@gmx.net true     true
3    local_admin  admins@email.earth    true     true

## Conclusion (and Other Platforms)

The Plain Linux approach works equally well on other platforms. I hope this helps.

Categories

## Numbers games: Z-Transform and some sympy

Back in uni me and control systems never really clicked. So I took a dive into the topic as I had the time. Long story short for all who know equally little or even less about the maths of control systems as I do:

• Often, you model control systems using the Laplace transform. That’s nice and good for continuous systems.
• Once you’re going to the discrete (digital) world, you want to use the Z-transform.
• There is a (complicated) relation between Laplace domain and Z domain.
• There are approximations, using a complex maps to map from Laplace to Z domain (a bilinear transformation, Tustin’s method, being the one used here).

This is a post on serendipity, because what this is about is actually neither control systems, Z-, or Laplace transforms but a numbers game. Bare with me.

Many systems can be modelled using a rational function of the form

You can easily put this into sympy and play around with it.

import sympy as sy

s, z, Ts = sy.symbols("s, z, T_s")

a0 = sy.S.One
a1, a2 = sy.symbols("alpha_1, alpha_2")
b0, b1, b2 = sy.symbols("beta_0, beta_1, beta_2")

G = (b0 * s**2 + b1 * s + b2) / (a0 * s**2 + a1 * s + a2)
G # use print(sy.latex(G)) to get the LaTeX representation

Nice. Let’s generalise this a little bit:

import sympy as sy

s, z, Ts = sy.symbols("s, z, T_s")

N, M = 2, 2

a = sy.symbols(f"alpha_1:{N+1}")
a = (sy.S.One, *a)
b = sy.symbols(f"beta_0:{M+1}")

G = sy.Add(*[b[i] * s**i for i in range(M+1)]) / sy.Add(*[a[i] * s**i for i in range(N+1)])

By playing with the values for N and M we can now generate any order of rational function we wish. Again, not the thing this article is about.

Then there is the bilinear transform called Tustin’s method that approximately maps the function G from the Laplace domain into the Z domain.

Let’s do this in sympy:

tustin = 2/Ts * (z-1)/(z+1)

G.subs({s: tustin}).simplify()

Sweet. Now, I want to have the representation in two polynomials in z (numerator and denominator) to extract the coefficients (because that’s what I wanted to hack into code):

num, den = G.as_numer_denom()
num = sy.Poly(num, z)
den = sy.Poly(den, z)

num.coeffs()
den.coeffs()

This gives you a nice, codeable result:

num: [  T_s**2*beta_0 + 2*T_s*beta_1  + 4*beta_2,
2*T_s**2*beta_0                 - 8*beta_2,
T_s**2*beta_0 - 2*T_s*beta_1  + 4*beta_2]
den: [  T_s**2        + 2*T_s*alpha_1 + 4*alpha_2,
2*T_s**2                        - 8*alpha_2,
T_s**2        - 2*T_s*alpha_1 + 4*alpha_2]

… but my eye spotted a pattern here. Do you see it? It’s easier to see in matrix representation, and I’ll only look at num (i.e. the expression with βi in it, because I set ɑ0 to 1):

Here is the resulting matrix for N = M = 3 and for N = M = 5:

Here are the resulting matrices for N = M = 4:

Here are the resulting matrices for N = M = 6

The odd values for N and M produce matrices with all non-zero values, so I’ll look at those first.

What I find stunning is the fact that there are so many patterns and yet I fail to think of a pattern that creates the whole matrix. For example, the first column clearly follows the binomial distribution:

where k is the row number starting with 0.

The first and last rows of the matrix are power series of 2 and -2 respectively:

where l is the column number starting at 0.

Every column has a different sign pattern starting at all positive on column 0 and alternating, starting with + on the last column.

Ignoring the signs, the columns are related like this

Maybe this follows from the power series mentioned earlier. But looking at the coefficients in each column the binomial pattern is lost as you move to the inside.

And last but not least wild zeros appear when you look at the even values for N and M:

Apparently, in the middle row as well as the middle column, every other value is zero.

All right, I’m absolutely sure mathematicians and control systems theoreticians have looked at this back in the day and this is all solved and done. Still, I found it interesting to look at.

If you want to create those values yourself, you can use this little script, and modify N and M. But beware, at N = 10 the script runs for 10 seconds.

import sympy as sy

s, z, Ts = sy.symbols("s, z, T_s")

N, M = 3, 3  # <-- Modify here

a = sy.symbols(f"alpha_1:{N+1}")
a = (sy.S.One, *a)
b = sy.symbols(f"beta_0:{M+1}")

G = sy.Add(*[b[i] * s**i for i in range(M+1)]) / sy.Add(*[a[i] * s**i for i in range(N+1)])

tustin = 2/Ts * (z-1)/(z+1)
G = G.subs({s: tustin}).simplify()

num, den = G.as_numer_denom()
num = sy.Poly(num, z)
den = sy.Poly(den, z)

for l in num.coeffs():
p = sy.Poly(l, Ts)
for i, B in enumerate(p.all_coeffs()):
try:
coeff, _ = B.as_two_terms()
except:
coeff = 0 if B is sy.S.Zero else 1
print(f"{repr(coeff):6s} & ", end="")
print()

Oh, and of course you can change N and M independently from each other, but I could not be bothered to look at this as well.

Categories

## Rogue docker processes on Ubuntu: Snap!

I’ve come across a situation where out of three docker containers that should be running, only one was shown by docker container ls -a. Weird, hu?

The missing two were clearly running their processes because the web service they provided was active and reachable.

After a lot of forth and back, it turned out that Ubuntu 20.04 with snap came with a docker snap – and had a recent Docker CE installed via the apt package manager. Don’t do that.

Once the snap docker was removed, it worked and showed all the containers again.

Categories

## Missing important details: Don’t update to DSM 7 yet*

*) if you rely on the SynoCommunity package rdiff-backup to do backups.

That’s really it: The package is not supported (yet?) and backups should continue.

Don’t repeat my mistakes.

Categories

## Repairing Electronics Rant

So, I have this Samsung NP900X3A laptop that I used for about 4 years. It started to have random glitches when I was about to start some very important work, so I decided to get a new laptop. I never sold or scrapped it, and recently I took it out of its neoprene bag, cleared the disk, installed a clean Ubuntu 20.04 and it worked like a charm…

…ish. Well, obviously it’s old, the battery had seen better days, but for work when it’s on power supply, it is a very decent device. Until recently:

I plugged it in, opened the lid, pressed power aaaand: it won’t boot.

Okay, my first though: The SSD died. Nothing catastrophic, a 120 GB mSATA costs about 45 €. So I booted into a Linux live system, quickly installed smartmontools (smartctl), opened my favourite smartctl reminder wiki page (seriously, man pages are still the best, but they’re so long!), and found that all tests passed. Well, the drive was old, yes, but still intact. fsck did not complain either.

But it still wouldn’t boot.

So, I opened the BIOS menu. First suspect: time and date were off by years. Aha! The BIOS battery, or actually, the battery for the real time clock (RTC), had died. And that set everything to defaults. Including: UEFI [DISABLED]. Well, you could’ve printed something, you nob…

Anyways: Enabled it, it boots again.

So, what about the battery? It’s obviously a coin cell on two wires with a plug. Nice idea, very accessible. Really. But at least some documentation would have been nice. No label, no print, nothing.

Measuring reveals: diameter 16 mm, thickness maybe a little below 2 mm. Hard to tell exactly while still in the rubber protection. So, could it just be a CR1616 or CR1620 cell with some spot welded soldering pins enclosed in some heat-shrink tubing?

No way to get the original part anymore. Samsung stopped selling laptops for quite a while in Germany (or even the whole EU, though it may have recently restarted again). Ebay has some offers. Anything from 5 to 150 € (seriously, it’s a coin cell, who charges more than 10 € for that?).

I decided to get one full package (cell, heat-shrink, wires and connector) for 5 € and a bare cell with spot welded solder pads for 2 €.

So, where’s the rant? Here:

COMPLEX ELECTRONIC DEVICES “BREAK” BECAUSE OF < 5 € COMPONENTS AND YOU NEED A F****ING PhD IN BATTERIES AND EBAY TO FIND REPLACEMENT PARTS BECAUSE NOONE THOUGHT OF PRINTING “CHECK UEFI SETTINGS” OR OF PRINTED LABELS ON A BATTERY OR A SIMPLE REMARK IN THE USER MANUAL (“YEAH, AND BTW, THE BATTERY IS A CR1616 FOR 84 CENTS WITH A MOLEX CONNECTOR, YOU’RE WELCOME.”)???

By now, this laptop would’ve been thrown out twice by every single average capable consumer in the world. And honestly, I cannot blame them. Right to repair will have a loooooong way to go.

So yeah: Turns out the coin cell actually was a CR1616.

Cut off the wires, placed heat-shrink, soldered to replacement cell (not to the cell directly but to the solder tails), pack all in heat-shrink and voilà: a 1,65 € (+some solder and heat-shrink) replacement for a 10 year old laptop.

Categories

## A Toy Traffic Signal with an Arduino Nano

My kid likes traffic lights. So I decided to quickly build one. I started off on a breadboard with a bunch of components:

• 3 × LEDs (red, amber, green)
• 3 × 470 Ω Resistors (only had two, so the third are actually 2 1 kΩ in parallel)
• 1 × 10 kΩ Resistor
• 1 × Pushbutton Switch
• 1 × Arduino Nano (actually a clone)

The circuitry is relatively simple. All LEDs’ cathodes point to ground (GND) through a 470 Ω resistor each to limit current and not blow anything. The anodes are connected individually to the digital port D3 through D5.

The button works such that, if open, it connects digital port D2 via a 10 kΩ resistor to the 3.3 V output of the Nano. When it closes, it’s supposed to pull D2 to GND, so it’s directly connected to it. The resistor serves to limit the current and – again – not blow anything.

That’s it. Now, for example, pulling D3 high (i.e. applying 3.3 V) will switch on the red LED. Likewise for D4 on the yellow and D5 on the green LED. Let’s introduce some human readable names instead of D2 through D5:

#include <Arduino.h>

#define BUTTON      D2
#define LED_RED     D3
#define LED_YELLOW  D4
#define LED_GREEN   D5

In the setup() function, I set the button pin as an input and the LED pins as outputs:

void setup() {
pinMode(LED_RED, OUTPUT);
pinMode(LED_YELLOW, OUTPUT);
pinMode(LED_GREEN, OUTPUT);
pinMode(BUTTON, INPUT);
}

To control the traffic light logic, I scribbled a Moore machine on a piece of paper, starting with a pedestrian traffic light (only red and green).

The machine starts in the red-light state. It remains there, unless the button is pressed. Then, it instantly switches to green (a pedestrian’s dream) where it remains until the button is pressed again.

This behaviour can be implemented by putting a simple switch case statement in the main loop() and implementing a state variable:

// States
enum states{ ST_R,
ST_G
}

void setState(bool red, bool green) {
digitalWrite(LED_RED, red ? HIGH : LOW);
digitalWrite(LED_GREEN, green ? HIGH : LOW);
}

void setRed() { setState(true, false); }
void setGreen() { setState(false, true); }

void loop() {
switch(machineState) {
case ST_R:
// Serial.write("State: Red");
setRed();
break;
case ST_G:
// Serial.write("State: Green");
setGreen();
break;
}

Why use this enum thingy, you say? You could just use an int and call your states 0, 1, 2, etc., you say? Yes, you could also sign up for the biannual global brainfuck contest.

But, hey, wait, there are no transitions!

Exactly! For that purpose I’ll use an interrupt.

void ISR_ButtonToggle() {
Serial.write("BUUUUTTON!\n");
if(machineState == ST_R) {
machineState = ST_G;
} else if (machineState == ST_G) {
machineState = ST_R;
}
}

And don’t forget to actually bind the interrupt function in setup():

attachInterrupt(digitalPinToInterrupt(BUTTON), ISR_ButtonToggle, FALLING);

Now that’s a simple and pretty useless traffic light (well, I still like that I can switch it to green at will, this really should be the default pedestrian crossing light!).

But it is easy to extend. You want a yellow state? Add a yellow state ST_Y to your state machine (the switch case statement). Or you want it to show yellow before switching from green to red? Add a yellow-before-red state ST_Y_R. Or you really like to let people wait? Add a wait state. Oh, and of course you can call your states whatever you want.

Here is my full implementation of a non-pedestrian red-yellow-green traffic light as it is stated in the German traffic rules:

• From red to green, with wait periods, go through
• Red
• Red+Yellow (yes, red+yellow, yellow meaning “get ready” and red meaning “not too fast, it’s still red, Freundchen!”)
• Green
• From green to red, with wait periods, go through
• Green
• Yellow (yes, only yellow. The meaning is disputed. Reasonable people say it means “Prepare to stop, if you can’t make it in time, it’s okay.” while idiots say it means “STEP ON THE GAS!”)
• Red (This is also sometimes disputed, I won’t go into this.)
#include <Arduino.h>

#define BAUDRATE      9600
#define BUTTON        DD2
#define LED_RED       DD3
#define LED_YELLOW    DD4
#define LED_GREEN     DD5

// Transition delays
#define DELAY_R_RY        250
#define DELAY_RY_G        250
#define DELAY_G_Y         250
#define DELAY_Y_R         250
#define DELAY_R_G         250
#define DELAY_G_R         250

// States
enum states{ ST_R,
ST_RY,
ST_G,
ST_Y,
WAIT_R_RY,
WAIT_RY_G,
WAIT_G_Y,
WAIT_Y_R,
WAIT_G_R,
WAIT_R_G };

states machineState = ST_R;
bool builtin_led = false;

void setState(bool red, bool yellow, bool green) {
digitalWrite(LED_RED, red ? HIGH : LOW);
digitalWrite(LED_YELLOW, yellow ? HIGH : LOW);
digitalWrite(LED_GREEN, green ? HIGH : LOW);
}

void setRed() { setState(true, false, false); }
void setRedYellow() { setState(true, true, false); }
void setYellow() { setState(false, true, false); }
void setGreen() { setState(false, false, true); }

delay(delay_ms);
digitalWrite(LED_BUILTIN, builtin_led);
Serial.write(".");
builtin_led = !builtin_led;
}
Serial.write("\n");
}

void ISR_ButtonToggle() {
Serial.write("BUUUUTTON!\n");
if(machineState == ST_R) {
machineState = WAIT_R_RY;
} else if (machineState == ST_G) {
machineState = WAIT_G_Y;
}
}

void setup() {
pinMode(LED_BUILTIN, OUTPUT);
pinMode(LED_RED, OUTPUT);
pinMode(LED_YELLOW, OUTPUT);
pinMode(LED_GREEN, OUTPUT);
digitalWrite(LED_BUILTIN, LOW);
digitalWrite(LED_RED, LOW);
digitalWrite(LED_YELLOW, LOW);
digitalWrite(LED_GREEN, LOW);
pinMode(BUTTON, INPUT);

attachInterrupt(digitalPinToInterrupt(BUTTON), ISR_ButtonToggle, FALLING);

Serial.begin(BAUDRATE);
}

void loop() {
switch(machineState) {
case ST_R:
// Serial.write("State: Red");
setRed();
break;
case ST_G:
// Serial.write("State: Green");
setGreen();
break;
case WAIT_R_RY:
Serial.write("Wait: Red->RedYellow\n");
machineState = ST_RY;
break;
case ST_RY:
setRedYellow();
machineState = WAIT_RY_G;
break;
case WAIT_RY_G:
Serial.write("Wait: RedYellow->Green\n");
machineState = ST_G;
break;
case WAIT_G_Y:
Serial.write("Wait: Green->Yellow\n");
machineState = ST_Y;
break;
case ST_Y:
setYellow();
machineState = WAIT_Y_R;
break;
case WAIT_Y_R:
Serial.write("Wait: Yellow->Red\n");
machineState = ST_R;
break;
case WAIT_R_G:
Serial.write("Wait: Red->Green\n");
machineState = ST_G;
break;
case WAIT_G_R:
Serial.write("Wait: Green->Red");
machineState = ST_R;
break;
default:
Serial.write("Woops...\n");
machineState = ST_R;
}
}

Is this the best way to write this? Hell, no! Is it even a good tutorial? You tell me. I think it’s fairly simple and shows some evenly common concepts:

• Using IO pins
• Wiring LEDs without burning them or the output pin
• Wiring a button*
• Using #define statements to make code more readable and lives more easy
• Using some abstraction (setState(), setRed(), etc.)
• Using a state machine as an extensible concept
• Using an interrupt (Why use an interrupt? Because it interrupts your code exactly when the button is pressed. It makes your little circuit very responsive and – at least in this example – it is even easier to write.)

*) You may encounter a problem that goes by the name “bouncing” (or in the beautiful German language “prellen”). That’s when your button press toggles multiple interrupts in a row. Try putting a capacitor in parallel with the resistor, or as a bad software hack, add a 5 ms delay (delay(5);) to the end of your interrupt function.

Anyways: My girl loves it. Ha! You thought she’d be a boy because it was a traffic light. Shame on you!