Performant, scientific computation in Python and Rust

The Rust programming language gained a lot of attention over the last years, and began to slowly infiltrate the Python ecosystem with an ever-increasing number of tools and libraries in the Python ecosystem such as Ruff and Polars which are implemented in this language. Unlike Python, Rust is a system language optimized for performance and memory safety, and some consider it the spiritual successor of C++. Despite its steep learning curve, it is the perfect candidate for extending Python and its ecosystem when performance matters, in a modern and memory-safe language.

This session demonstrates the path of creating a scientific package in python (following best practices and modern tools) and gradually migrating parts of it to Rust for additional performance gains. The use case is a naive implementation of the "Expectation maximization for Gaussian Mixture Models" algorithm from scratch, a relatively simple yet efficient machine learning method. The session addresses the following points: How to build a Python package with a modern tools set, how to translate a numerical algorithm into vectorized Python, and optimize the package with a performant Rust implementation of the critical parts. Prior knowledge of Rust or the algorithm is not required. Note that the goal is not to learn Rust in this single session (this requires at least three days) but rather to provide a superficial overview on what makes this language so great and well-suited for extending Python.

Participants are advised to follow the clone the repository below and follow the installation instructions to avoid longer download times during the session. https://github.com/StefanUlbrich/PyCon2024

This session took place in track Programming & Software Engineering and was classified suitable for intermediate domain / intermediate python by the speaker.

Transcript (auto)

Auto-generated from the recording utilizing Open-Source AI. Speaker labels (Speaker 1, Speaker 2) reflect diarization, not identity. Timestamps refer to the recording.

Speaker 1 [00:06]

Yeah, thank you. Welcome everybody. Sorry, I've got a bit of cold. So if I cough a bit and it's too loud, apologies for that in advance. Yes, first of all, who has ever written anything in Rust yet? Okay, that is quite a lot, actually. Who uses a tool written in Rust? Yeah, that's expected. Probably there's people not aware of using tools written in Rust, especially Python community. Yes, I'm very passionate about Rust. I started about two years ago writing on private projects and now can do it in my company, too, which is really nice. And I'm also trying to convince a lot of people of going into Rust and Rust is cool and we have a small group and we're gaining a critical mass now and, yeah, that compares to a cult. which is why they made a better title slide for me, with the help of artificial intelligence, and I gave this talk internally earlier, it was called the rust in the Python, it reminded me so much of a 90s point and click adventure game, but I had to keep that slide, apologies for that one. So yes, first of all, who am I? I've been introduced, thank you very much, I'm a data scientist at Blue Yonder, we have our booth up there. Feel free to visit us later. And so my job as a data scientist, I take probabilistic forecasts in retail and compute the optimal order amounts to make them produce less waste. So it's a pretty fulfilling job that I have currently in, and I can write it in Rust, which is like a huge bonus on top. My research interests, I still do research with university, like on a hobby side. I'm a roboticist and interested in in cognition and machine learning, especially self-organisation, and that's a bit related to the topic that we will have also today. So I started Python way over 15 years ago, so I consider myself experienced on that field. Rust is like two and a half years now, and the last half year I do it as a full-time job. So I hope that qualifies me a bit to give you some of the experience I made with those two languages. So what to expect. The idea of this tutorial session, and for me that's the very first time I give this format, so bear with me, and it's also the first time that people queued up for listening to me, which is really new. So it is a time lapse of a person or a development process where you're given like a research paper and you want to put it in the end into a Python package that is highly optimised even in Rust. And all that process we want to squeeze in in this one and a half hours, I hope that we succeed with that, so that is the main idea of that, that you get like a whole view of it and for, I mean, the people who are not coming maybe from academia and haven't done that and might be a bit intimidated when they see papers with a lot of formulas, I want them to give them a bit of, well, confidence that not, be not intimidated by just a bit of few formulas, does not, should not scare you off, yeah, it is doable, it's all scientific numerical or machine learning package, you can do that, right? So we look at how this process will look like, and then you can decide for your own, and hopefully you agree with me in the end. We'll look at the strengths and weaknesses of Python and Rust in the process, I probably don't go so explicit into that in the slides, but ask also, yeah, at any time feel free to interrupt me and ask questions, right? And how well those two languages play together. Opinionated, I try to give my best practices in coding style and packaging, but there's as many opinions as there is packaging, right? All is open source. On the slide before, I did mention there was a QR code going to the GitHub repository of this talk. It's also been in the program. Everything there is open source, including these slides, which are in pure HTML, so you You can watch them in real time, basically, and scroll back if you found something. You can use also the whole repository as a template if you like. It's open. Equally important is what should you not expect. This is not a comprehensive Rust course. I linked to a comprehensive Rust course of the same title from Google, which in my opinion is the best way of learning it these days. And I gave that in the company also, so it's like a five-day course, and that's probably the time you should at least calculate to get a certain level of confidence. Might be less, because if you just, like we do here, take some part of Python, make it faster and put it back, then you maybe don't need the whole picture. But to learn the language, you should have five days, especially if you have a background at C++, that would help. No deep learning frameworks, so we will work with pure NumPy on the Python side and with equivalent on the Rust side. Likewise, no complex models. So there won't be LLM, there won't be transformers, there won't be neural networks. So I chose something that is simple to understand, that I can explain, and yes, just to keep it in the timeframe also. There will be bugs. This is a live session, the first of this type. I might try live coding. I tried that once and felt badly, so maybe I chicken out. Everything is prepared, so I've got everything in its own branch. I can just go through the code if need be, and if you find typos, keep them. So to the structure. Structures, as I said, the first part, and hopefully can't keep that short, will be the boring theory part. So it will present one special algorithm that is very dear to me, which I have probably coded for like 100 times in less programming languages than that, I guess, but quite often. I know it by, yeah, I know it very well. And I think it's well suited for given its complexity and so on. So on the second one, we'll switch over to Python. So we have now these formulas and the description of the algorithm. How do we bring that in efficient Python, that would be the second part. I try to do the timing as good as I can, not spending too much of these first two parts, because the last one will then go into Rust, where I'll give a brief overview, so people experiencing the language might be a bit bored in that first part, it's just a repetition, of course. We'll write basically one part of the algorithm in Rust, and package this as an extension, which then will be called from Python. And a tiny bit of benchmarking in the end to wrap it up. And that should hopefully conclude the tutorial session. Again, feel free to interrupt me if you have any questions. So for now would be a good time, especially between those blocks, I'd say. Okay. So let's first start off with the Gaussian mixtures and expectation maximization. Just maybe a quick poll. Who knows that algorithm? learning and is familiar with that, I guess, okay. K-means? Better. Perfect. Gaussian mixtures and expectation maximisation. K-means is learned by expectation maximisation, so that is not new, and Gaussian mixtures are simply the soft version of it, literally speaking. We'll get to that in a second. Why this model? It's easy to understand. I think it's educational but still yet useful, so you can do real-world problems with that, and I used it in different contexts actually for productive environments. It's fast, it's versatile. We implemented in less than 50 lines of code. I think that's also helpful, right, if you can grasp it. And it's equally easy to train and evaluate. So I don't want to get a large language model corpus of words and wait five minutes until it's trained, I'll just generate random blobs and that's it. Bonus, it's probabilistic. I think that's very important in machine learning because a probabilistic allows you explainability, so you can justify why you have results. Why not? Yes, as I found out doing that, creating that, well, tutorial, also when I prepared it the first time, I realised it might not have been the best choice because you can express that algorithm fairly well in array notation, or in matrix notation, which means you can write it in NumPy, PyTorch, and so on, which, again, are in C, and, yeah, fast. We can see there will be an improvement, but it will not be drastic. Imagine if you have something that involves caching, a dynamic programming thing or something with a non-linear of workflow, control flow, then you will get speedups that are drastic, 100 times or more. I mean, you see it from if anybody uses UV, the new packaging system, which also adapted for this, they are like easily 100 times faster in solving the SAT problem of resolving dependencies, right? And that is because it's written in a language that supports that kind of speed. There is a question about the slides, have they been put in the repo? I can hopefully demonstrate it in case I still know, I lost my internet connection. Sorry to interrupt. No, fine, let's see if I can connect again to my phone, if not it is hopefully still open. No, it's not. That's probably even the wrong... Great. If I had internet, which I still have, apparently, I don't know why I do have. Come on. You see it on the right. You just click on that link, and you'll end up in the same slide. It's just the GitHub pages. It's all sponsored by GitHub. If you're interested, it's Reveal.js. I opted for trying out something new and not doing PowerPoint. You lose a lot of time if you do. But it's kind of cool. I've got a lot of formulas in it, and that's important. The formal definition. Let's rush a bit. A Gaussian mixture is the distribution of a points, where each point is drawn from a single Gaussian. So, imagine you have these kind of clusters, it's useful clustering, I can spoil that, and it's a probability distribution that first decides is it one of these clusters, and if so, it draws from that Gaussian. And that is often enough to do a clustering in not so high dimensions usually. We get to that also, too. So the equations are pretty straightforward. As I said, it's one of those, so it's like an or. Or in probabilistic, it's a sum, so you've got a sum that's weighted by a probability that the model is the one you're going to pick from, and then it is just normal distribution, and I just put the equations there, but we won't even see them and Python code. All right. There's an alternative formal definition that is, if you put it in so-called plate notation, that would be more like probabilistic modelling, people would do that, and, yeah, important is we operate on data that is partially observable, so we have got, like, 2D points that we had before, that would be the X, and we have got a set, It's a variable that we can't observe, so it is from which Gaussian did we do that sample. We don't know that. If you just see grey points, they were nicely coloured. If you don't see that, it's difficult to tell from which it is. Also every time you train it will have a different number, so you can't tell that, and it's basically called a latent variable. That's also why it is unsupervised learning, basically. We have one categorical distribution, which is a big dice, and you roll the dice to get out which gaussians you take from, and then you have a fixed number of normal distributions. That's like in k-means, so you have a fixed number of clusters, but it's a probabilistic model. You can expand it. If you put the Dirichlet process on top of it and make it really complex, then you can and scale it up to infinite numbers. Examples where you could use that would be take an animal distribution where you look at height and weight. Those will be correlated, usually, and that's dogs, and you would have different dog breeds and clusters on that. Or if you have a demand forecast from the domain I'm from, you will simply sell more steaks in the barbecue season than you would do in winter, right? You would have different distributions. That is a combination of both. I wrote a small blog post. It's LinkedIn here. So I've got a small blog which I neglect for at least a year, but there's at least a few blog posts on it, and I wrote about that one without formulas. Okay. Gilson, you train with expectation maximisation. Expectation maximisation is an algorithm mainly for clustering and unsupervised learning, and exactly for the case where we have missing data, we don't know in our case which goes and produces the data point, so what are we going to do? We call that repetition, we call that latent store, and we have got an alternating pattern. We have got two steps. Then an expectation step, which is basically we guess that value. We start with an initial assumption and we guess what the set is, and based on that guess we improve our model. Then we guess again, then we improve, then we guess again, and it's alternating until you converge, which you will, but it's unfortunately also pretty prone to local minima, but that's one of the problems that you have with machine learning, of course. You can train more than that. You can use that for k-means, you can use it for a certain type of neural networks, or self-organizing networks, which is kind of my speciality, other mixture models, hidden Markov, or Kalman filters. So it's pretty versatile. And I created a small animation with the code that we're going to create today. It's running a bit too fast, I guess. But what you can see, which distinguishes it from K-means. So K-means would look more or less similar, but all the circles would be, well, they would be circles, would be circles, not ellipses like they are here. And each point would have the same color. Here you can see basically each point has a blend of a color. That means when you do the expectation, you compute for every point the probability that it belongs to one of those possible clusters. So it's called a soft assignment. That's what I meant, it's the soft version of k-means. Small question. This one is chosen, well, I chose for the example, I think, 10 clusters, and per chance I picked the right number of clusters for learning. Yeah, the question is, of course, legit. Yeah, you never know how many clusters to pick in front. You would usually pick more than that. When you train, you will see that a few of them become irrelevant, or you can go with Bayesian criteria and so on. Personally, I had that problem, too. I went into the full de-resolving route, which simply means you assume that you have an infinite amount of them, and that it's just highly unlikely that you have a high number that is relevant, and then you would converge until the correct version. But that's definitely a bit out of scope to code here. I think I said everything about that slide already, so we start with the model here indicated with the two means and For each of the points we give it a color which would be in that case Yeah, two values that sum up to one because they're probabilities and they indicate what it is in the end We end up with a matrix called the responsibility matrix because we have n points and we have k probability Possibilities where it could belong to so we have got an n times k matrix In formulas, this does not look so intimidating, I guess. You again have the probability of belonging to a model, so for each column of that vector you would do that. You have the probability of what you would create when it came from that model. That's just Gaussian, and then with Bayes you can get a probability. So before that you have got densities, they are not normalised between 0 and 1, but this is also a normalisation step, so you get a real probability in the end, and the equation given, so it is really not that intimidating. A point of, if you want to use that, that red minus one is the bad thing, because computing an inverse matrix in high dimensions is in cubic complexity, and that's why if you look up in sklearn, they say it is not suited for high dimensions, which is not completely true, because it means if all your variables are allowed to correlate, so you can trick there also in going in lower dimensions, so if you just have a diagonal matrix, then it is actually pretty fast even. The maximisation step is way easier, I would say. Simply it is you have those weights or basically those probabilities computed in the step before, You know how much a point belongs to which cluster. And then you just have to compute a new mean, and a new covariance, and a new weight. So three values you have to compute. And they are actually straightforward, because there is just the sample mean, but weighted. So formulas are equally simple. So where you would simply have the sum of the x divided by m, you replace it with the weight. So you get a weighted sum, and divide by the total amount that you would expect. The same is for the covariance. That looks pretty similar, actually, the equation. And to get the probability that you're in a cluster that's like simple counting, more or less. It's straightforward. So this algorithm is pretty powerful, given how simple it is, I think, and also how simple it is to implement in the end. So that is, again, the justification for the choice of that. You can use it in higher dimensions. Again, the warning, it will explode at a certain point. But there are ways around it. Also, the higher the dimension is, well, the larger your space get. And that means picking right starting values becomes non-trivial. But there's also means for that. But it becomes definitely more a challenge, like everything becomes. And you can also use it for regression, which is actually pretty interesting, that you can use an unsupervised machine learning method for supervised learning. I won't go into that one, but it's actually pretty cool, because you... I will go a bit. Sorry, I like that topic. If you have a complex function, you can basically, when you do that, partition the surface or the latent space that you want to explore into compartments that are mostly linear. And it will make that in an optimal fashion. And then it's easy to do a regression step in the end. The regression is simply in the inference. So that can help if you don't know how to partition your input space and so on. You can use that for linear models. So it is powerful. You can do quite interesting tricks with it and with little effort in the end. But these little efforts mean you have to code it on your own. Circle again. Why are we doing this here? Okay, the Python implementation. We have now spent nearly half an hour, so I'm a bit below speed. We can still try at least a bit, start with an online coding. I will first, yeah, a few words first. I hope you're all prepared if you want to follow the steps on your own machines. I invite you, if you don't want to do that right now, I completely understand that, But feel free to use the github repository as a playground, right? So every significant step is its own Its own branch so you can iterate through everything and see whether you did something wrong or not And it guides you basically through it. So that was the idea of the whole of the whole scope and yeah Giving you a bit a guardrail of doing that new task task. Design decisions. I kept it simple. It's not like something I would put on PyPy for sure. I've got a small data class which has a numpy array for means, which is of the shape. It has the number of components and the number of dimensions. It means each row Represent one of the bubbles or centers of the of the of the gaussians Coviances is pretty clear. That's the times D. Which is cut off at the end. So it is K times The metrics are it's good to have a convention there because your compiler even in rust will not help you I'm telling you if you mix up the dimensions and that's already always the problem when you do it in live coding also We'll have a function that does expectation to be filled out, so it takes the learned model and it will take the data. What it doesn't do, it's not a method. I want to be able to freely exchange it just for the sake of testing. In reality, that would be a self and it would be a method, of course. But for now, it is not the main focus of it. It's more functional programming in that sense. We have got a maximization step, still not implemented in the first step, of course, And we do a very, very simple initialization process. It's actually one of the simplest, but the best and the highest dimensionality, in my opinion. That means simply we compute for each point a random responsibility. So we just roll it, not a dice, but we decide how much each point belongs to each cluster, but per chance. It's just drawing from a Dirichlet distribution, basically. Benchmarking will just be done with timeit, but I put that at the end in the end. So the expectation step is what we begin with. It's probably not intimidating. I think I went into the equations already. What's important is in Python, if you use a lot of control flow extensively, so if you loop. If you use a list comprehension, as cool as they are, they are loops. That means every time you operate in the Python world, you enter in the C world of NumPy, you go back and everything you do there is going through the whole overhead of taking care of which type it is, checking whether methods are there, and everything that the interpreter that takes away from you, which is nice, actually. Python is cool because you can just start coding. You don't have to think about types. And if you're thorough, you just add the type hints afterwards. Probably you should start with them. Yeah, it gives you a lot of freedom. And yes, I would say for prototyping, personally, I'd be at least five times faster in Python, at least, compared to fighting with the borrower checker. The type system in Rust is really, really pedantic. But if we can avoid loops, we should. So we try to do it without loops. And the only interesting part of doing that is in computing the covariance matrix. So we have this equation. You see we do have already a few loops in there. So basically, you have to do both the denominator denominator separately, and the sum in there is basically the loop, and you've got a bunch of loops in the Scala product hidden as well. So we can try to first put it into a matrix notation. We can get rid of the sums and just imagine how we would have it laid out in memory. That's something if you're doing Rust, you're always thinking how stuff is laid out in memory. That's addictive. And you can put it basically in a form where you've got a matrix with three dimensional arrays, the one is transposed in the first two dimensions, and the one in the middle is just as a factor in between. And that at the bottom should pop up later, so that's the fun of doing it. So the letter, basically this notation is called Einstein notation, so Einstein apparently had to fight a lot with indexes and his physics stuff, and kind of invented a way of writing that down, which basically, already in Python more or less, would look like that. You have something with that three dimensions, something with that three dimensions, with that one, they all have the same k, I want the k still be there, but I have the n which which disappears, so I want to reduce along that dimension, and the D still remains there. And that's enough, actually. Yeah, that works. You can describe what you want to do there with the Einstein notation. And how cool would it be if we had an Einstein sum notation function in Python where we just put in a string and put in the parameters and it works, and, well, we're still speaking about Python, aren't we? Of course, we do have. It's a numpy. It's easy. It is there. This one won't work, because it does not really like. It's difficult to explain. Showing how it will like a work. So sorry, too fast. Oh, I didn't even fix it here. OK, it will not work, because those two d's are not allowed to have the same name, because they remain here. And C would not be able to figure out which of the same Ds would be the first and which would be the second. So you have to replace one of those with an E, which I forgot to put on the slides, but these are the things you'll find out very quickly because it complains when you do it. Yes, I think I skipped the live coding part, if that's okay for folks. That's embarrassing for me for now at least. We go through the maximization step first, I'd say. We can, maybe if there's a bit of time, we can go through the expectation because it's way easier, actually. Okay, what we do is first checking whether shapes are all right. Then we will need the sum of the responsibilities along of the access to the zero. Actually, let me look in my cheat sheet whether it's correct. That means you recompute for each cluster all the responsibilities. right? We need that to find out how relevant it is, right? If none of the points is claimed by that one cluster, it will have a straight zero here. So this is basically the denominator that we will see quite often, I think. For the means, it gets a bit more interesting. So if we are working in that high-dimensional thing, we have to kind of blow up objects in order to be compatible. So if you look at the data, the data has n elements times the dimension. Why the responsibilities have n elements times the number of clusters. So these two are not compatible. You could not multiply them. So I probably repeat it for most of the people here, but you know that that's broadcasting. In Python, you can simply say, please, NumPy, do a copy in that direction. It actually doesn't do a copy. It's smart enough to just pick what it is supposed to do. But if you went to, I guess it's in PyTorch, it's the same. I'm not completely sure. I used QPy, which is like a drop-in replacement for NumPy. It would literally copy those dimensions, blow up the thing and waste all the memory of your graphics card. But then it would be tremendously fast. On your PC, you have to do it in a loop. So what you do is you tell it by adding a new axis that this one is to be repeated in that direction. And you do that on the right spot. Then you will end up with an object that is compatible. And that gives us then, in the end, the means, which are of the shape k times dimension. So that's what will come out because we reduce along the first axis, which makes we get independent of the number of points. Of course, we don't want our model to grow with the number of points we put in. So that is the obvious choice to remove it. It took me at least a while to wrap my hand around that way of thinking. It's easier for you, hopefully. But, yes, that is the way how you would write efficient numerical stuff in Python. Covariance, we just had it. As I mentioned, here again I substituted the later D, although D and E are equal, for that algorithm it would not need to be, at least not for the formula, so that's why you have to assign two different values to it. But it doesn't change anything at all. In the end, the weight is simply, yeah, counting how much weight or how much influence we have in total for each of the clusters. And that's how you would basically do it in Python. So no live coding, I guess, but we can still quickly jump to, there's the final code, walk through everything at least, be thorough, if I find it should be an init, I guess, yes. So initialization, as I mentioned, is we want to assign for each point that we have, basically. For each, yeah, exactly. I create basically a data set and no, no, no, no, no. I construct manually a responsibility matrix and I call maximize on top of that. what I did here. Just to have a valid... Oh, there's a question from the mind. Is it more possible now? Great. One more? Ooh. That is huge. Happy? Okay. Cool. Yeah, So basically, I create a responsibility matrix and call maximize on that. That will guarantee that I have all the correct shapes and initializations that I need. It is trivial. It's simply to have the thing running. The expectation is probably the more interesting part, and it's not optimized at all, but I did not run to write down the Gaussian equations again and again and again, which is way too often in my life, I guess. I just pick, basically, the multivariate normal that is in SciPy, call the PDF on it with the data, just iterate in a loop over my... It's comprehension. Thanks. Ah, perfect. That's the ideal audience. With a list comprehension. I was looking for a word. Put it into an array, it's not efficient at all, it's not meant to be efficient at that point either, just to make the algorithm running. As I did that once in life, and I found out that you make mistakes which you can't detect. Quick question in the audience again, who does test-driven development? So who writes their test really before they start coding? You really do that. My deepest respect. Cheers for that. I couldn't actually imagine how you would do it any different. I started doing it. It took me more discipline than I probably have. But it helped a lot. But that was when I came into a new domain where I here at least I think I do have some expertise and there I didn't have. But it would have worked way better if I had in the last session had some unit tests because they would immediately have told me, no, you can't do that. These areas don't fit. And so you have at that stage, the moment you check out the unit test branch, you will get unit tests that give you at least a small insight. They They are probably not very good. They do their job, especially the test maximize, where you basically compare the values given a trivial data set. The covariance metric is where I compare them element-wise. And you see the tolerance that I have is huge. But it is enough, at least, to fit for the live coding, which I think was enough for that stage. I guess that would be it for the Python implementation, if I'm not mistaken. Is there any questions at that part? Yes, there's one. We do an initialization when we call it a maximize function in the end of the initialization, correct? Yes, I hope. Maybe, yes, here we call it a maximize function. No, we don't need it. It's just for initialization. I do a dummy expect by creating manually the responsibility. The question was, why don't I call expect? I just wanted to have valid objects, to just have a starting point. I did not really care too much about that, to be honest. If you did it right, so if you're really interested in how to implement that correctly, You would not have an initialization at all. You would have in the very first expectation step what you do, you would start with that. So you would do that in the learn loop kind of, checking whether it's the first iteration. If it's the first iteration, you would do the random assignment and then continue with the algorithm. Here it's really, the focus is not on that function. The focus is simply bootstrap the algorithm and let it run, which it does.

Speaker 2 [37:12]

It's better because it's recorded.

Speaker 1 [37:17]

So it's actually just an example in Python house algorithm can be Written, but it's not working one because we are going to write the working one in rust. It is working I would say it's working. I would not say it's beautiful Okay, but when maybe can you one more? I still don't understand where I wear the expect function called There's no expect function at this and the since I bootstrap with a random initialization Yeah, and the random initialization is not in yeah, but you okay. Let me Start over What you see in in blog posts and well in in textbooks you would pick random values for the mean You will pick random values for the covariance and for the weights with a high probability that a few of those are near where close your samples, especially if you are in very high-dimensional places, and I usually am in high-dimensional places. What you can do is simply, I say these points, they get already an emulated expectation step, and they have already the values attached to it. Then I call maximize and maximize will read the all the cluster will be huge and You remember maybe the animation at the beginning started with huge circles that went small. That's exactly because of that and The advantage of is you cover all of your data. You will not neglect any of your data point No one will fall off the table. So actually I'm thinking about that. It's probably the better way of doing it. It's not a cheating way of getting it done. It's be where you should probably do it even. Let's take it offline if there's still questions on that. If that's okay for you. Thank you. That's interesting. I love having discussion about that. All right. We have quite a bit of time left. That session goes until half past. Okay, then we could have done live coding, I guess. Maybe we'll finish a bit earlier even, I don't know. Let's try. So, okay, now lay back a bit, it's the Rust introduction part, it's addressed to the people that are new now, yeah, so who haven't seen it and just heard, oh, this magical language and I want to, everybody talks about it. First, what I skipped, this is Ferris, so that's the mascot, we have got the nice snake and they have got a nice crab. And that's also the reason why we had this crab-cold image, if you have not figured that out, of course. It has been developed around, well, start development started at Mozilla at 2010, approximately, And the first stable version arrived at 2015. So I remember reading about that, but I don't want to learn a new hip language, but then we had a very strong voice in our company, and that convinced me to try it out later. What it definitely does, and I agree on that one, it feels modern. It feels like really modern. It's not like if you go back to C++, you have got that feeling of the 90s, even if you haven't lived in the 90s. You will feel like in the 90s immediately. Which is probably not fair. C++ went through a lot of cycles and iterations and became way better and also copied some or got inspired by some of the concepts in Rust. But as Rust is from 2010, it started way later. And the later you come, you don't have to repeat the mistakes of the others. So they could throw a lot away which was not useful, and C++ has to stick with it because it has to be backwards compatible. A very good example is actually if you call a function and you don't write anything in front of the variable, then this variable can't be changed. It is not writable. You can only read it. In order to change it, you have to write short for mutable in front of it. So you have to be explicit about whether you want to change it or not. Actually there is a small bug in the code before, I can't quickly show that, actually. I didn't fix it because it doesn't really matter at that stage. But I take here the input data, and I blow it up, oh, I didn't talk about that line earlier, and I subtract the means and make it bigger, and that's basically the base I use. But I think I modified the data here, don't I? I'm not really sure. Do I not? That's shattering? Okay. Then it's not a bug. Yay. But it could happen easily. I do modify the GMM, definitely, yes. Which is what I intend. But it's not explicitly stated. Yeah, and that is sometimes the reason for very, very weird bugs. Yeah? You smile. You've seen that. But having it explicit is way, way, way, way, way more better. So that's just an example of it drifting away again. It's a very active community. It's very committed. I ran often into some things that I could not explain to myself. Also because I tend to pick problems that I want to put in a library and to make them as generic as possible and really pushing at least my limited understanding of the typing system. And you write somewhere in Discord and you have an answer five minutes later. So as equally amazing as our community is in Python, I would say. Smaller probably, I'm assuming, but it's continuously growing. Most of you probably know Hacker News. Just search for anything but written in Rust. That was like a joke for a while, yeah? And people complaining, oh, you did something new, but you have written it in Rust, yeah? But it's showing how the language is growing, yeah? People are really adapting it, want to do things in it. Rewrite it, and there's an overlap with Python, yeah? I mean, look how many tools are written now in Rust for Python, especially in the tool chain where really waiting is annoying, yeah? Where Rust is amazing. I sometimes have to enter mistakes in my code just to see whether it's still alive. It's so fast. Polars, last year we had at least three Polars talks on this conference, if I remember. UVE is pretty new now. I used it in my project and it was amazingly fast. Where I would wait 20 minutes, it does it in a second. Also showing the fact that there's so many tools written for Python shows that the communities are also aligned and overlapping. There's an intersection or joint set which is pretty large. For me, and I linked there to a blog post which inspired me a lot, it impacted also how I code in Python. So I'm using the type system in MyPy very extensively, and I would always use protocols, unions and what so on. And the impact it had on my Python code, I mean it comes from the C++ world, I'm used to typing, and I think Python is a bit dangerous because it lets you slip through so many things that makes it complex. If the code base grows and you work together with other people, then that quickly can become a problem if you're not disciplined, and I'm for sure not disciplined. I know that. Yeah, but it influenced me, so I really adapted patterns, and I had one pattern also in here in the Python code, but I'm drifting away again, so probably not going there. But I would recommend reading the blog post, it was really like a small eye-opener. In Rust, we have object-oriented programming, but we don't have any inheritance. And for sure, we don't have multiple inheritance, and that's a very good thing. There's huge discussions whether inheritance in general is a good thing or not, and I don't want to get into that flame war, there is simply none. What we have instead is trait interfaces. They are the closest thing that you would see in Python in protocols. I don't know, has somebody used protocols in Python for MyPy? Good. Very good. Proud. They are amazing. And I think it helps from the concept. So a protocol or trait or in Java, it's called interface, is basically describing the behaviour of something. So it does not tell you it is inherited from something or inherits something of it, it just describes the behaviour. So you can basically make an interface, you say anything that implements for instance summary has this method, and you can separate it from the implementation. So you can have that in a completely different file, I mean if you remember the C++ days you would have header files and distribute basically of what the person sees from the implementation, which I think is sometimes very good, because I looked at the GMM code from sklearn, and that one is bigger than 50 lines of code, like, way bigger, and part of it is in Fortran, and I don't know, it's huge. So having that separation is nice, right, and you would then basically implement for that trade summary for the news article, so that would basically be the syntax. You will run into things that are a bit difficult to digest at the beginning, just because it's now on the slides, I'll just say it. You can't have functions with an arbitrary number of parameters. Right? So, like, the print function can take as many, you've got the ellipses, the dots, or the arc, star, and whatever, you have the possibility to put anything you want into Python, of course, and it won't complain. It will complain at runtime. But Rust needs to know how much memory it has to reserve on the stack. So on one data structure that grows, the more you go into functions. So it needs to know what comes in, right? So you don't have the opportunity of putting arbitrary amounts of arguments into it, which I just did here. And that one you go already into metaprogramming, which is in Rust like something obscure, which you would not do. In Rust, you do have that. They are called macros. And you will recognize them always with an exclamation mark. You have that in the print. You have that in the format. You have it in a few convenience functions. We'll also get to that. But, yeah, for now, it's simply that's the way how you do a function that expects more than a fixed number of parameters. That sounds limiting. I did not miss it, actually. Rarely does it come to that cases. And if we had that, we would actually write a macro. That would go too deep, but there is a simple way and a very, very complex way for extending the language, actually. You can write plug-ins for the compiler that lets you do basically anything you want. So it has still a lot of freedom review, but it's getting more and more complex. Speaking of complexity, one of the things that everybody talks about Rustin, which is probably the one biggest selling point because even the White House put out an announcement a couple of weeks ago, is that it's memory safe. In C++ or in C, you can create an array or a memory of sequence with five elements, and ask for what's written in the 100th element, which doesn't exist. You just go into a different memory and look what's there, and that's used as entry points for cyber attacks, for instance, right? So I guess there's a number of 75%, a large percentage, which is actually due to misuse of memory or overflows in languages like C or C++, I have to say. You can avoid that, you have shared pointer, you don't have a garbage collection, in C++ you would do a reference counting which would avoid that, but still it's very, very easy to mess up and to create unpredictable code. What Rust guarantees you, you will never get out of of bounds. You can get out of bounds, but then the program will panic. So the program will fail. Hard. But you know where. You can't violate memory. And there's a different point. You can't do unpredictable behaviour. That's simply not possible. You can't have a race condition in threads, in a threaded program. That's going too deep. I'm personally not very good at multi-threading myself, so Quick example here the the tool for that is is it's called the borrower checker. So We have some sort of fixed rules that you have to comply and if you follow these rules there will be guaranteed that these things can never happen and Why the compiler enforces these it's usually It's not test driven development. It's compiler error of driven development in rust It tells you exactly where to look, it's fun, try it out, it really tells you, nah, that doesn't work, try it like that, and usually it's the right answer. It works really, really nicely, right? So quick example, we would create a mutable variable. We create a string, we get into a context. So that's basically a limited of the lifetime. So yes, we use curly brackets. We don't have indentation and rust. Think of that as a good or bad thing. It's better for anonymous functions. That's for sure. If we leave that, our run would go out of scope. But we put all the value of I have to check it again. Yes, and the variable still works. So in the end, probably I didn't bring that over very well, you have to take care that if you have got a variable or piece of memory, that there can be multiple places or threads that can read it, and no one can write, or one can write and nobody can read it. So you have to think about how long does something live, and you should not be surprised if things go away, and you have to do something against it. You will get always a good note. The good thing is you won't have implicit copying. In many languages, you would just make a copy and you wouldn't even notice that. In Rust, you have to be very explicit. If you want a copy, either you give it a carte blanche and that's implementing the copy trait, or you have to call explicitly clone. And that means, yeah, imagine you have got an algorithm that calls one loop millions of times per second. If in that time you make a copy of something all the time and copy the memory move it around get You think about a lot about how the cache is organized. Yeah, you have got this huge memory. I've got a small cache in Close to your CPU, which will make things very fast You won't don't want to break that yeah, so rust helps you go being closer to the metal I would still say it's fairly abstract, but you have to think about these things or you can at least and And really squeezed the last bit of performance out of it. I got a bit drifted away on that slide apologies So enums enums are great. They're amazing. They're much more powerful than they are Basically in every other language, I don't know whether somebody has an simple Similar concept. It's they are actually algebraic data types The closest thing and it's a bit like that is if you think of them as union types in Python So, you could have a function returning a union and that would more or less fit what you have there. You just don't have only a value. You can have values, but you can also have if messages of type move, then it has that data attached to it of type XY. If it's of type write, or instance write, then it would have a string attached to it. If it's color, it would have a triple of integer values. And in memory, it's organized that it will have just the largest necessary memory block reserved to represent all of those. And all error handling, all return values, whether it's like the optional in Python, are done with that system. There's no exceptions in Rust. You don't have an interruption of the control flow of the program, which is much too expensive in terms of computation. We always return an object that tells us everything is OK, here is your result, or it is an error of that type, and that uses this one. Once you get used to the concept, it really, really, really starts, you start to appreciate that. Much more. You have got pattern matching, aligns very well with enums and what we have seen before, destructuring. You have got a lot of functional programming, hopefully we get to that in a second. And the tooling is probably the best what you can get right now. I mean, if they speak about tools for packaging in Python and dependency management, you always will hear Cargo as a benchmark, right? And that's the tool for us. So it's really, really nice. And all those things, the pattern matching, everything, the functional programming comes at zero cost. So that means while it makes your programming very comfortable, it is not slower in any way as if you did for loops or in any way else in a classical C-style programming, it would be exactly as fast because they don't put an overhead. I am so baked into the language that that works. Okay. I guess you noticed that I like that language. I know. Yeah, so now we're going to port that MSTEP into Rust. I will stick to the slides a bit for that first. Just to see how that looks like. It's not Java. You don't have to create an object just to have a main function, luckily. That reminds a lot of C or other functions. You need one main entry point if you have a program that's always a function, main, brackets are mentioned already. If you have a variable, or you want to create a variable, you have to put a let in front of it. You would put let mute if it's mutable, and here, just also to show that it is operating well with NumPy. I just read data that I've got in one of the notebooks in the last branches. You can create a small training data set for the unit tests, and things can fail. In Rust, if something can fail, you have to react to it. I mean, you can call something in Python, and if it throws an exception, well, it's up to you to put a try and error block around it, right? It's your choice of doing that, and you can easily forget it. In Rust, not, because this thing does return an error, and in order to get to the value, you have to tell it to give me the thing if it worked. If you're lazy, you can simply say unwrap, which means give me the value or die, so the program will crash if something happened. Fine in unit tests, everywhere else, probably not. Aero handling I won't touch today, so it's unwraps everywhere, but you shouldn't do that. Again, print is like the Python statement, a new line, you can't put arbitrary stuff into it. it has to be a macro. I've said that before. Rinse and repeat for the other two things. So that's a bit looking familiar, hopefully. Yes. If you want to pass data into a function, it's getting more complex. Even the simple things are a bit more complex. I agree on that one. There's way too many colons in there. I think that is... It should be two. Sorry. I don't know why it's four. So we have here a function that takes an array of dimension two and returns it. I think that code either way would not work right now, but it's not intended. So the way it's written before means basically if you do that, the caller that puts in the data into that function will not use it ever again because it gave it away. you give away the ownership away and you don't make a copy, you just, if you wanted to use it again later, you get an error. So you moved it away. If you want to pass it in without losing the ownership of it, so here you have it, but give it back later, yeah, then you have to add an ampersand and that is just like you would do it on the C word. It's called a reference. What you do in Python, like always, where you pass values always by reference except for basic data types. And that's the same in Rust. Basic types are copied. You have to do that that way, basically, to tell him that you don't want to copy it. Now that got a bit away. If you want to be mutable, I've said that before, you have to a few typos in there. I'm sorry. But then you have to put the before. And this example basically takes the ownership of a view object, which is fine, which is only temporary. That one you would not have to put as a copy because it's just a view on something. That's a little overhead of copying. And assigns it to the content of the other. Put the values into the one. This one would basically overwrite the result. that's why it has to be mutable, of course. Array passing, that's about returning. That will fail, definitely, because you have always to give back a data type, and you'll see quickly what I meant with why Rust is a bit complex when it comes to the typing system. So in In Python, you do have generic data types and generics in the type system of the whole language is generic, but you do have that. In Rust, you use it extensively. You would never write a function just for one data type if it works for all, then you would make it generic. In that sense, I take here an array which can be either float or integer or string or or whatever, whatever is implemented for being able in an array, but it has to be dimension two, that one I want to have be fixed, and my return type is kind of pointed to, quick check, so I never point afterwards. Yeah, basically, I fix it to float at that point. Yeah, okay. I quickly skip it at that point. Error handling is probably the more interesting point because now I compare now Python to Rust and how basically if you use the right libraries, there's a library called ndarray which mimics the interface of NumPy very closely. And yeah, I go through a few examples examples how similar it looks like, right? So if I compute the responsibilities in Python, I would just sum the responsibilities along x is 1. Yeah, it should be both 0 or 1. If I do that in Rust, then it looks pretty similar, right? But I have to put in an x as object. But you see the familiarity. You have got the let in front of it because you are doing a variable assignment, which is mandatory. And interesting is broadcasting. I mean, that's what we do quite a lot in our code that we looked at. Yes, it looks more complex. I agree on that one. Problem is also if you try to address, yeah, a slice of an array, you are kind of accepting that you have arbitrary parameters. right so that's why we see a macro again but they use it so often that they give it only one letter to keep it short trade-offs I guess and yeah if I want to do though the calculus that I had before I would call basically on the two broadcast items the multiplication operator but everything is a function again so it's the multiplication is a function that you implement somewhere And if you call that you will lose the memory of that data Yeah, so you put your objects into the data and they're gone for you You can't work with them anymore unless you put an ampersand in front of it and signal I want to use it continue using it, right? so That's very explicit agree on that one, right? It's You can make that wrong. You will make that wrong at the beginning. Sometimes it will be easy to resolve But the thing is you can't not cannot accidentally copy things. So you're kind of forced into the right lane of Implementing efficiently Dot product basically the same we don't have the at operator So we have to call transpose dot that's on again. I don't want to lose X. I want to use X further around the route, and I want to use the parameter down the route, that's why the ampersands are there. There's no Einstein sum, unfortunately. That's really sad, unfortunately. Which would make the code a bit more complex. Great. We have, like, 15 minutes more. That's stuff a bit. Prepared, of course. I did not really talk about the structure of the whole package which I would love to do quickly now. I've got even a cookie cutter for that. So if you want to really start from scratch, it's also my GitHub repo or just copy this This one, of course. You have all the code of Python in Python directory. You have all the code for the Rust implementation in Rust. And you have a different folder for the bindings in the basically the connection to the Python world and the bindings folder. Everything is defined in, yeah, that's pre-commit, of course, but everything is in PyProject in there. So a bit of tooling in good style. I'm using Poe the poet as a task runner also. That is maybe some tricks to look into. Maybe copy something if you want. And UV as the dependency resolution. That's all written in the readme that hopefully works as it's written there, but that sometimes does not. I hope it does for most of you at least. Great, that was the excurs. Going back to the Rust package. Rust has always a cargo.toml at the root. We can quickly open it. It just tells the package name, versions, dependencies, and that looks, I think, very familiar to what you would expect in Rust, right? At least the PyProject.toml is new enough, and they probably drew from each other inspirations, I would assume. You have also something for examples. I have data, of course, there, and I've got a benchmark, which I use later to show the difference of speeds which given the time you're probably not able to touch at that point so I created here and that's basically when you go to benchmark or if you go to maximize one, two, three, four I've got like four versions of how you could do the same thing they are working basically starting with one where you quickly go over them where you would iterate so that's how a for loop looks like. Hopefully you don't see them too often in Rust because loop sucks. If you can do functional programming, you would opt in for the functional programming all the time. So that's how a for loop or a range looks like. You see that it is kind of optimized to be modern and not too complex to write. I remember the C++ days where you've got to its.begin until its smaller end and E++ and whatever. That does not exist, right? So this one is fully functional code. It resembles a bit of what we had before, but it basically does the Einstein sum in a loop, and it works, but we can do it nicer, I guess. Definitely what you have to do before, and that's something that you also would have to, would be better in Python, that you first create the array with zeros, just to reserve the memory, because every time you want to blow something up once it has been already created, having a list and adding elements to it is ridiculously expensive. You have the arrays, you can do append and Python, you will do that all the time, but it's really ridiculously expensive to, oh, I'm running out of memory, I fetch somewhere a bigger memory and I copy everything over. If you start working in Rust, you pay attention to that. You don't do that anymore at one point because you are there in order to be efficient. You wouldn't be there if it were different. I quickly jump to the last one, which is a bit, let's see where I go. So there the interesting part is I've got the covariance created and adjusted. I didn't call it like that, but that is basically the, you take the mean values and subtract it from the data. So this cleaning of the... I don't remember the term in English right now. So you subtract that, and you have got the responsibilities. I take all those three, and I want to iterate over the axis, one in the first case, in the second case I want to iterate over the second axis, and I want to iterate over the second access in the responsibilities. Here's a mod because I want to change it, because that's my target where I write into, right? It has to be usable. That's already great. And it won't compile, but let's just hide that for a second. I go into that, I get basically a tuple over which I can iterate. I get an object that is an iterator which has in it items of type tuple or triple in that case. For each of those, I will do, and that's what is a lambda in Python, sorry, I don't remember the name already. It's called a closure. I think I've heard the term also in the Python world. Closures capture what's around them, like they do in Python. In a closure, you can refer one of the variables and it will be put, it will be given to that object, and that's quite interesting when it comes to ownership and what happens, that's way too deep right now, but it is powerful, especially compared to C++ where you really have to say I want that, I want that, I want that, and you have to write it down because it has been added later. Here, and I do here the regular thing they would have done with the Einstein sum. This code, I think it's not that complex for now, I agree, but usually you would see this pattern like you have an iterator, you would do a map, you would do a filter, but a filter is nice, but you would also take while, you would just cut the elements at the beginning. Well, skip while, you would cut the elements at the end, and so on. It would look a lot like what you do in Pandas basically, like if you have got this chaining of commands, and that's the way you would write in Rust in the end, and then this becomes so readable, because it's simply, okay, this step, this step, like a pipeline, you can map so many real-world problems into it, and it really becomes a charm. It's lovable if that's at that stage. And this translates to exactly the same machine code or program code as the first example of the loop. Those examples are trivial. Of course, it's just one loop, and here it's just one iteration, but no matter how complex you get, it will be the same code as above. So you won't have any difference in runtime. Sometimes slight, because reality is a bit different always, but neglectable, I would say. So that's cool, yeah? And iterator, it's functional programming, it's nice, and what is really cool is if I put in that line, because if I have that, That way of writing something in a functional code opens the door for parallelism. If I have that line, it will immediately run on all of my code, of course. It will split the workload on all of my calls of mine. It doesn't have that much. It's pretty old. But it would put it at least on the four cores that we have. Immediately with one additional line of code. And that would give you a speedup factor in that example. Not too much, but it's at least factor three that we have here. Just one additional line. That was, for me, at the first moment, pretty impressive, actually. Yes, question? Yeah, a question of understanding. Is it right that this is possible for just the reason that there is no data sharing? So that you really give away the responsibility or the ownership of data? Yes, there is a bit more to it, of course. So, the inter-par iter is not something of the standard library in Rust. It's a package called rayon, which is pretty popular because it is so freaky easy to use. But numpy has to cooperate, and not numpy, sorry, and the array has to cooperate. So, these objects, the arrays support it, that they can be split up into different parts and iterated each independently over it. Yes, because the area is one memory block held by one variable and then this would not be possible, right? Yes Exactly. That's not a complex workflow. Yeah It is parallelizable, of course But it's pretty cool, and what's also cool, and I definitely won't do that in that session, but we could look at the assembly code now, and it's again pretty old, but in newer machines you would see that it uses AVX, SSA2, and what so on kernel extensions, right? You can do that explicitly, and unfortunately that's not in stable Rust yet, but it will hopefully come soon, that the portable SIMD comes. So SIMD stands for Single Instruction Multiple Data, so that if you can do one thing with one float, you can do it with five most regularly, and you would have a wider register. That's why Intel and AMD gets more complex every time, because they get more registers for that, instead of a sensible vector architecture. And, yeah, so that's amazing. You really get all the nice benefits of being close to the metal, and still abstract a bit from it, right? It's still, you start thinking in where's my memory and how do I squeeze something out? Okay, you, that's why you're here in the first place. But it still gives you the comfort of a very modern feeling and language. Okay, we are approaching the end, I guess. I've got a few notebooks with it, so you can, where did it go? You can basically, if you want to, the tutorial at home with more time, of course, You can run through it. Most of the code should be functional, I hope. If you find a bug, don't keep it. Please give me a pull request or an issue or something that I can fix it. But quickly run just functional test and just to see the speed up, basically between the Python version, because that's the thing why we start the whole thing, right? So it's benchmarking, the last part, actually, of this talk. The speed up of having this pure Python implementation, which should be fairly optimal because it does not do anything in Python that much, right? The code was very, very... It still has to read your instructions. It doesn't build up a graph like PyTorch would. Yeah, PyTorch would be faster. No discussion about that. That would be more optimal. And I doubt I would be able to beat PyTorch code here. But NumPy has to iterate through everything I've written, so it's a bit slower. And you should be expecting some impact at least. Not that much, but a little bit, and to cut it short, so we would have like 13 seconds with a pretty wide, no, not pretty wide variety, but a factor 10 at least variance around it. And the Rusco that we have seen, well, is a fourth, only 25% of that speed. So that's substantial, I'd say, not that much. As I mentioned, I would recommend do that on Levenstein distance, For instance just to warm up something that is dynamic programming where you really have some More complex things tree operations. Yeah iterate over graphs or trees. That's something that you can't put in blocks easily. Yeah, then It will be much much much more I guess in a speed up nevertheless. I think it was a Entertaining example at least hopefully Bottom line for me, Python is high level. It is very quick in rapid prototyping. It is very good if you can put it on the shoulder of something that is very optimised. So machine learning, data science, I mean, we tried running Rust in a notebook. We won't try it again. It's not fun. Python, Rust is not out there to replace Python. That is something I can at least guarantee and it won't. It's simply a different thing. It is not, you would not, you don't have to write everything in Python because in a few things it is simply slower. But in the most real world applications where you use Python, it doesn't really matter. And the important part you just offload it into that domain, right? So we do have five minutes, don't we? Ten. That's easy. Okay. I skipped half of my talk. Not really. Not benches. There will be a dedicated session about PyO3. I think tomorrow. I don't know. It's right now. Then I probably can just hand it over to the next one and make this one short. I, of course, use PyO3. It is the defective standard to communicate between the world. I completely skip the word unsafe. Rust forces you to operate in a way that your code will always be predictable and will not violate anything, no memory, no boundaries, and nothing wrong can happen. Unless you tell it, no, ignore what I'm doing. That is what you call unsafe. We can put a block of unsafe which allows you to do some things that you're not allowed to do usually. not completely a carte blanche, but it allows you to talk to something else, for instance. If Python wants to give me a memory, I access memory that doesn't belong to me. I'm supposed to do that. So I say, yeah, let me do that, please. And PyO3 does that for you. So usually I've never written unsafe code. We did it in the team, I guess, but I personally haven't. I'm not looking forward to it. You would have a lot of boilerplate code here, as to expect it. I think in Boost and the older times from C++, it was more, I think, as I remember it. But it's fairly straightforward. You create a Python module, you've got a lot of imports. I define a function. It's a bit tricky with picking names. It doesn't find it if it's not the same name as it has in Python, I think. I tell it it's a function, the function name, I accept Python object which represents my compiler or my interpreter, I expect two data, and I return a triple of data. And from these, I just convert them so I can use them, so basically you get something from Python, it has to be extracted what I want, that is done by NDArray for me, and I just call my function, and I'm done. We don't even have a return statement, right? If just the last thing and it doesn't have a semicolon, it's implicitly the return. And I return everything went right. That's it. Still bigger. But I'm done. Okay. Yeah. That would be the... I don't go deeper, especially for the next session, the very next session. Benchmark I mentioned a bit. Slightly different values here. 30 milliseconds. I did not prove it because I did not run it, but I ran it quickly before that session. That's it. So key takeaways. Be confident in implementing something scientific. Especially maybe for the younger audience. Just because it has a few equations shouldn't scare you off, right? It's nice. Just pick one of the algorithms you like and just implement them, right? Even the complex Gaussian processor or something, yeah, complex formula, it is ten lines mostly. It's rarely more. And it's fun. Yeah, and it's fun also creating mixed Rust Python packages. I hope I could convey that a bit. I really like that. It's cool. It doesn't always make sense, and you should also be aware of that. don't optimise out of what is not worth optimisation. There's another key takeaway from me, joke at the end, never ask an AI to generate a picture with a crab and a python together, or you get that one. It doesn't spit out a single good picture. There's always horror movies like that. There's a question from behind. Sorry, you had your hand up quite a while. Apologies for that. Just a quick question on the benchmarks. I don't know if it applied everywhere but I saw that, I mean the NumPy memory layout is basically role based and I saw that quite a number of column based operations. Did you ever try the other layout and see how Python behaves then in the benchmark? It would definitely affect the caching, like a lot. So yeah it would be, I did not go for full optimisation here. To be fair, it is like a year old, I guess, the code, when I started it. And then that was still a hobby for me, right? So it was really spare time, and I did my steps with Rust, I still consider myself a learner after two and a half years. I'm not an expert by far, intermediate maybe, with a bit of luck. But yes, you're right, different member layout would be better, not only for caching, also for SIMD. That actually makes a lot of impact if stuff lies next to each other. So you can move everything in a single instruction. That makes a lot of difference. That toy example, you probably wouldn't have seen that. You would need a bigger number and crunching, but there would be a substantial difference, yes.

Speaker 2 [83:47]

We still have five minutes and there is one question from Slido. If you want to put it on Slido, please put it. If you want to raise your hand, we still have time. So how to get started with our own first Rust package? Any recommended boilerplate?

Speaker 1 [84:09]

Well, I've got a cookie cutter, it's really difficult to pronounce that word, that you can use as a starting point. You can use that GitHub repository as you please. That's what it's there for. It's MIT. You do whatever you want with it, of course. The interesting part is more, you know, you've got no solutions, so you have to look for good problem. That's more the relevant part. The question is what to do with it. It's not like I would be able to get the idea, oh, no, I'm writing a SAT solver and I do the next UV thing. Well, A, it has been done, and B, I would not even have thought about that. You have to find something that inspires you to do something. For me, it is picking some interesting algorithm, mostly related to my research that I still do on the side. And, yeah, that's how I did it. I can only speak for myself, of course. But I think the scientific approach, that's why you're probably in the talk, would pick something that is simple. Pick something easy. Pick self-organizing maps is great. Gaussian process is great. Something really bottom.

Speaker 2 [85:23]

Or maybe Einstein's submission convention. That would be interesting.

Speaker 1 [85:28]

Interesting, I've seen somebody trying starting with it, but it ended What if I want to multiply a matrix on GPUs for example, I don't know if it's possible weak point Good point on weak point that you found. That's a bit of a pain point GPU You let's be honest you speak about in video that isn't C++ there's a little way around. I personally looked into a bit of that and that would be the route of WGPU, so WebGPU, which gives you the possibility to run on in shader, compute shaders. It won't be as efficient as it would be with Nvidia I guess but it's at least a start but you wouldn't mix languages because shaders right now would have to be written in W, GL, SL, in that one shader language for that, and there's a project that tries to establish Rust as a shader language. That is very promising. It's by a game studio somewhere in Europe, I think. I don't remember where it's from. It looks very promising. I'm watching the repository. I'm always hoping there is something coming up. There was at least 10 approaches. Rust is still in that phase where everything pops up. The common joke is there's two games written in Rust, but five game engines. It's evolving. You'll settle on one horse, it might die. That happens. In Python, it is a bit more established, but still can happen. Still things get obsolete, and the one developer does not want to do it anymore in their spare time. That happens obviously everywhere, but in Rust it's still a bit fresher, and GPU is for me still a pain point, yes.

Speaker 2 [87:23]

Okay. One more. What about debugging? Can you step from Python into Rust code?

Speaker 1 [87:30]

Deep learning debugging debugging. Yes, of course. Yeah. Yeah

Speaker 2 [87:34]

OK. Yeah, could you please also, I mean, I have a question from my end. It's really horrible, this creature, crab plus python. I'm sorry. So the covariance matrix that you shared was basically had some subscript as m. Was that like multivariate or univariate covariance matrix, like when you were just before the slide? Yeah, this one. So as I understand, this

Speaker 1 [88:06]

M is an index over the components. Sometimes I switch between that. I did not want to use integers. It just means one of that set of gaussians that you have. Model. I picked it for model.

Speaker 2 [88:24]

All right. Thanks, Stefan. Again, let's thank Stefan from all of us.

Stefan Ulbrich

About — in the speaker's own words

I am a researcher and programmer with a passion for geometrical methods especially folding and unfolding in machine learning and robotics. My interests are computational cognition, robotics, bioinformatics, machine learning and data science. I am an experienced Pythonista and Rustacean

Social card for talk: Performant, scientific computation in Python and Rust