Stefan Karpinski (Keynote): Julia for Data Analysis and Beyond
Thank you thank you all for being here and I'm grateful to be invited to. I keep getting invited to these talks at a PI data and Python conferences and I have to say I've always found the Python community of you really welcoming and warm and supportive and in terms of growth mindsets very much not threatened by the the success of others so I think I think we can say that python python is a growth mindset kind of community so i have sort of two positions i am two hats i am one of the cofounders of julia computing. Which is a company providing support and training and consulting for julia which is increasingly used and you know finance insurance aerospace all sorts of things. I also have a part-time position at the center for data science at nyu. So if anybody is at NYU and interested in talking about you know Julia using it getting involved let me know. I'm pretty easy to find on the internet so so I'm gonna start out with sort of a little bit of background motivation. What my data science stack looked like back in two thousand nine when I was doing data science but didn't actually yet realize that I was a data scientist. Because i think the term was kind of coined like right around then or at least. That's what google trends tells me so in a single project. At the time i was using matlab for numerical linear algebra and machine learning this was some singular value decomposition and non-negative matrix factorization stuff. I was doing and MATLAB is a really great tool for this so I that was pretty good and of course. I was in grad school so I don't have to pay for it. So that was okay too. I was using. RF was just a statistical analysis and for visualization because uh the visualization of matlab is just nowhere near as nice as our and of course the stats stuff is not not nearly as good either but then there were some things that needed to be written in c because they needed to be really really fast and efficient so i had some see in there and then of course i was dealing with lots and lots of data files lots of inputs.
Lots of outputs running simulations doing all sorts of stuff shuffle it around a cluster and that any of that stuff is a nightmare to do in any of these languages so I also had Ruby to tie it all together and I realized at this point that I had created a Rube Goldberg machine. If you're not familiar with Rube Goldberg. He was an artist. I think active in like the 50s and 60s and you draw these incredibly complicated machines. Where like you know this guy is eating and his spoon pulls a string which like flips a piece of toast up and then the parrot eats it and then that like makes the parrot go down and then it put like spills a bucket of water and yada yada yada and eventually like something useful happens. But you know it's just like cobbled together out of all of these crazy pieces and that's what i felt like i was doing now it wasn't using Python at the time had I been using Python or head Python been in the state that it is now back in two thousand nine maybe I would never start it on the Julia project. But but we'll see i'll talk a little bit about that in a bit so this is actually the for language problem but we sometimes talk about the to language problem and the to language problem it actually goes under another name that's older which is austere holtz dichotomy and this is after this was a term coined by John Ulster Holt. Who was actually the creator of tickle. Then he was familiar with that programming language and he observed at some point in like the 80s. I think that there seemed to be to like classes of languages. That were quite difficult are quite different and that were sort of. You know what things gravitated towards and they were systems languages and scripting languages systems languages are statically. Typed typically scripting languages. I Dana dynamically typed systems. Languages are compiled scripting. Languages are interpreted systems. Languages have user-defined types. Like you have a struct or something like that that you can define and then do stuff with scripting languages tend not to let you do that and instead they provide a couple of very rich standard data types that you can do all sorts of different things like link a dick right in Python.
It's this one like really really powerful data structure or pythons arrays which are actually the get called lists. I don't know why that seems like a mistake. That should've been corrected point but um you know you can use them as a cue or a DQ or stack or all sorts of like. It's a very powerful array data structure they also systems languages also tend to be fast whereas the scripting languages are slow systems. Languages are hard whereas scripting languages are easy so this is this is the dichotomy and this was sort of the state of affairs in like most most of like the 80s and 90s and now in the like to thousands. We're starting to see that dichotomy sort of come apart a little bit. We're seeing this this. The systems languages sort of take some of the the properties of the scripting languages and vice versa. So because of this dichotomy there's a standard arrangement that people have typically used to do things. Because you're like well you know. I want convenience. But i also want performance i want both of those things and the compromise is that for convenience. Use a scripting language at the high level. Something like python or our matlab but you do all of the hard stuff in the systems language so when you're using numpy for example like the the fast stuff the stuff that's actually doing stuff with arrays is all written in c or actually. Fortran right because you're calling gloss if you do a matrix multiply which is still written in Fortran. So that and that kind of works like this is a pretty good arrangement. This has been sort of the state of the art for 30 years or something like that. I have myself written many high-level wrappers to low-level functionality and then you know you're happy because you can you can get the performance but you get this nice ease of use to but it is not a perfect arrangement it's very pragmatic but it's got some issues so one issue is aren't the hard parts exactly where you could benefit from a better language right so like why are we doing the hard stuff in see that seemed like kind of perverse right like I actually did try to.
I was fixing us one of our benchmarks because I noticed that C was cheating and doing less work than all the other programming languages and took me like two hours to get this thing not to crash. I don't know what I was doing. I had forgotten how to do see index incorrectly and see and I was just like how god I'd forgotten how difficult this was. So so this is this is a real issue. This is you know. C programming is not super easy it also forces vectorization everywhere right. There's this sort of like you're trot pot. That like for loops are kind of immoral like it as it turns out there just like impractical for people it writing libraries and high-level languages. They're not actually a moral right. It's fine to write a for loop the only issue is it's going to be slow. If you're in one of these high-level languages you want to like push as much work as possible into the low-level language but sometimes it's really awkward and you end up writing some very strange convoluted code to try to do this. It also can often force like additional memory allocation right because you have to go through a lot of create a lot of temporary erase in this style whereas like in C or something you would just like use one array and then mutate it in place and that would be fine and you would not ever create any temporaries. There's another subtler issue that took me a long time to recognize which is that it creates a social barrier between users and developers right so in these two language systems the users use one language they use the high level language. The developers actually don't spend that much time programming in the eye level language.
They spend a lot of time programming the low-level language and use it being a user of the high level system doesn't make you better you know qualified in any way to work on the internals you kind of there's just a barrier there right and I even fight even though I'm like perfectly well according to my story not that good at C programming but like I've done a lot of c programming. I'm okay at it you know if I'm using Python or are I tend to like hit that sea wall and then I just stopped because it's too much of a pain to find out where the sea source code is and how things work and I just kind of even though I could. I just don't so maybe maybe that's the moral issue maybe I should but but people don't it's just not practical so Julie is a programming language tries to crush and shatter. This oyster holds dichotomy when. I we want to have our cake and eat it too. So you know these are the two columns and we do this. We sort of take mix-and-match pieces from both columns julia is dynamic. But it's compiled we let you have user-defined types and standard type. So the standard types are pretty powerful but you can also do your own user-defined types and I'll show you an example of that in a bit. It's fast and it's easy. How fast this is sort of the obligatory performance chart this is time on a bunch of sort of you know micro benchmarks. That aren't really super interesting. Except that they time things like how good are you. At recursion iteration scalar loops parsing integers printing stuff. You know doing a little bit of matrix multiplication. It's time relative to see the the y-axis is logarithmic. So you know being higher up is much worse and the languages across the bottom. Julia Fortran go JavaScript. Python Mathematica our matlab octave python is sort of where you like. That's as good as a good interpreted languages ought to be. That's sort of where you are but then you have these like compiled languages like Fortran and go and then julia is the odd man out because it's sort of its dynamic but it's actually within a factor of two of see consistently beaten people actually see the Sun real-life code to this isn't just contrived from Mike Reuben.
I think this is a much more interesting chart. This is the same y axis but the x-axis is actually how many lines of code normalized between whatever the minimum and the maximum one was for each benchmark each languages and you kind of see these two regimes. You have you have a long. You know you have this corner right here. These are the fast piled language. These are the system's languages. These are the scripting languages. So you have this kind of choice you like saw you can be up on this end of the curve or on this end of the curve. You can either be concise and slow or you can be fast and kind of verbose. Of course there's one point that's sort of in the corner here so that's where we want to be. I like I like being in that corner. So this is my data science stack today. Julia for numerical linear algebra statistical analysis. Yeah he's Julie for everything Python can give you a lot of this right like you could use Python for everything but the c line right so it gets you down to from four line which is 22 but to get all the way down to one language you need to do more and project sleep. I are sort of like pushing in that general direction so maybe in the future people will be using pi PI and like writing Python code and it'll be like fastest see so that's exciting. Um so for the rest of this talk. I am going to give you some examples. I'm going to do a lot of live. Coding and you know things might go horribly wrong hopefully they won't. I'm gonna demo how you can use julie as a general perfect language so not just a language for numerical computing. Although i'm going to give a numerical example but you know whatever you can also use it for data analysis the usual stuff like data frames and then plotting some stuff.
And then i'm going to show what hopefully is a compelling example of a superpower for a user and then Anna compelling example of how it gives superpowers to library writers too and so that's that's that I'm going to do in the Jupiter notebook so everybody familiar with the Jupiter notebook this is a Python audience so I assume so everybody knows about this right this used to be called the ipython notebook but now it supports like 40 languages or something like that as kernels. Julia I think was the first like real usable back a colonel aside from Python fernando pérez. Invited me and jeff and some of the other julia contribute contributors to come to berkeley and write this kernel and we went and did it was a fun week. Ok so this is this is what Julia code looks like you can do things you can say that max iterations is constant and then you define you can have a doc string it actually goes outside the function instead of inside that actually makes sense in. Julia but I won't get into why and you notice that I don't really put any types or anything on here. I just say you know the dude. This is not this is the Julia set that I'm going to compute. Doesn't really have anything to do with the name. But whatever it's a cute coincidence so you take two complex numbers. Z and see. And then what you do is you do this iteration step this Z you say Z equals Z squared plus C and. Keep doing that until it escapes. Beyond the radius of 2 and the number of tot like iterations. You have do before thing escapes. It gives you like a color which you then can visualize it and make pretty pictures which we're about to see so you just you know press shift enter and that evaluates the cell and you can see that you got a generic function now call. Julia with one method. And we can you know. Call it on one complex point and you can see. I'm a little too tall for this podium so I'm going to put this up here. We can see some other points like okay. That's a higher point. I don't know so takes more iterations to escape we can do this.
With a comprehension of two-dimensional array comprehension syntax. Here so this is similar to pythons one-dimensional array comprehension syntax but you can do arbitrarily many dimensions and get in this case a matrix so this is a matrix of the escape times across a certain range. I chose the range so that at least some of the numbers were interesting. And they weren't all twos. But you can you can see this is 501 by a thousand thousand and one alright. Well that's kind of a you know a lousy way to visualize things so let's actually load some packages that do more interesting visualization the first thing. I'm going to do is define a color map. Which is it's the color this. This package colors provides a couple of different built in color maps with this one's going to involve red and blue. I'm going to get a hundred points because max iterations was 100 so I'm going to go from 1 to 100 and then I'm going to use that color map to turn escape times. I see if you evaluate the color map. It prints it nicely. You can see the colors it's actually just an array of colors but the we know how to print that as a nice thing and you can see here now. I make an image out of it and you get this pretty picture. Can people see that reasonably well on the projector. I think okay. Cool one of the interesting things about the Julia set is that you change these numbers a little bit. You change the starting point a little bit and get kind of radically different things and you get a lot of different really interesting pretty pictures but you know going in there. Oh another interesting thing about it is that doing it that way is the Julia set if you take this. So this is essentially a four dimensional shape because it's too complex numbers as inputs which are each two dimensional. If we take the diagonal of this set it is the Mandelbrot set so that should probably be familiar to people so the Julia set and the Mandelbrot set are related in this kind of this kind of way but it's interesting to sort of explore this set so let's use this thing called interact which gives you.
You can take anything that so that expression. I I used to generate an image. I can actually just say put this macro in front of it at manipulate is a macro in Julia and give it a range for sliders and then we're actually going to be able to like you'll see in a second you can scrub the sliders around and get different images all right so now we've got that and we can change it. Which is this is a much nicer way to explore different things and you can see that it kind of gives you gives you different patterns and their / they're all quite beautiful so this is a little more sluggish than I would like so one of the things this is doing is this is generating a new array every single time and going through in you know and then constructing an image out of it this is sort of inefficient. It would be nice if I could just pre-allocate one chunk of memory and just redisplay that every single time and not have to like you know go through all this memory allocation because that's expensive so this is a version of that which is slightly more verbose but does exactly that so I use a let block to say that what these INR are so that I can actually use their length. That's just kind of nice so what I do is I create a data array and here you see one of the particularly Julian things. We don't shy away from talking about types so classically dynamic languages try not to talk about types. They sort of try to pretend they don't exist. We decided that you know the first thing. I mean the whole like the whole initial point of numb pile is to have like typed arrays in Python. Which is the Sun typed language and so in numerical computing inevitably. You end up having to talk about type so let's just have one good system for talking about types and make that universal so it's actually just baked into the language and the language knows about types so here we're constructing a two-dimensional array a matrix the element types are the elements are actually this.
RGB type with an unsigned fixed point. 8-bit value in there. But you don't have to know about but this is all sort of like custom types that are just defined in libraries but they actually have the exact representation you would like for the machine which is you can represent. RGB values is just like you know triples of aid but of 8 bits. Then you create this image and then you just sort of go over and instead of instead of constructing the thing every time we just take a for loop and go through and just update the image every time for loops are not evil in Julia. At least they're not going to be slow. They're actually going to be quite fast so now we can see that we can. We can scrub through this and it's significantly more responsive which is quite nice cool any questions about that so far all right cool all right. So fractals are cute but you know we're here for analyzing data not analyzing fractals so let's do a little bit of data analysis. I'm going to load this package called data frames which provides a data frame data type. You know very similar to the one that ships with our and very similar to the one that pandas provides so. I'm going to load this. Jared jared is here somewhere this is one of jared's data sets that i actually you know just played with at some point was like oh this is a nice data set. It's good for displaying them demonstrating things NYC housing it's got a bunch of columns about like different housing sales throughout the city. And you can see that the type of this thing is that at the data frame and then. I'm going to use this package. Called gadfly which is similar to the it. It's a grammar of graphics style visualization package in pure. Julia and when I say in pure Julia I mean like the whole thing is written in Julia there is no C code involved which is kind of nice. It means you can hack on it more easily but it from a user's perspective it kind of doesn't matter so ok.
The first obvious thing you want to do is. I'll get rid of the headers so that we can see a little more. I don't know plot some stuff so like square. Gross square foot versus estimated expense. That should have some sort of interesting relationship. Right when we look at this and it's kind of this like diagonal blob and you're like well yeah. I guess maybe that has a relationship so I don't know the first obvious thing to do here is like there are really big values and there are really small values and a lot of stuff concentrated in the little small very small value. So let's you know let's scale it to be log-log scale so this is this is how you add attributes. If you're familiar with ggplot this is sort of similar to the way ggplot works. Ok now. We're getting somewhere there you go it's you can. It's like pushed pulled everything out of that. Like far left corner. We'd we now actually have things kind of spread out in a nice clean way. Still not like a little bit would be nice two separate things out so one of the things you can do is you can color by one of the fields and I'm going to color by borough and you can see that now. Ok yeah color by borough this is getting a little interesting and we can zoom in a bit and you can see that yeah so. Manhattan is the blue up across the top and Manhattan is like significantly more expensive. You're sort of higher up for long this log log scale then you know say Brooklyn is next and then Queens and then the Bronx and Staten Island. I mean you're I bet a lot of people here from New York City this makes sense. I mean so you can do you can do other things. Pretty straightforward stuff like plotting a histogram of you know so since we we figured out that like you know its income and square footage thing or have this log log relationship you know if we take a relay to take a ratio that should actually give us give us some like fairly consistent. Stuff so we can take a histogram plotting by color of these things and you can see again like you know there's there's kind of two modes for Manhattan and then Brooklyn is sort of further down branches.
You know further down is this the same same data just sliced a different way. So yeah you can do data analysis pretty straightforwardly. I I think that's you know that's useful. It's obviously interesting and like necessary that you can do this but this is not one of the things that sets it apart because you can do this kind of data analysis in any of these tools so. I'm going to talk about some of the other things the things I think. Give you real like super powers so the next thing. I'm going to talk about this is super powers for users this example this guy. Patrick Usul dinger. I assume that as a German name he emailed because and we get this a lot this happens is pretty typical like I wrote Python code and I poured it to Julia and it wasn't faster. What am i doing wrong or you are a liar because you said this would be fast. Depends i prefer the you know the former but the other thing happens to and unfortunately our mailing lists are full of people who are like love hacking performance right so this is always like this is always kind of fun because everyone is just like. Oh yeah here are several like cool ways you can make this much faster and like frequently. It's like orders of magnitude faster by the time you get done so in this case it was an interesting journey and performance because they're a couple of like the original so the original. Python version took point five seconds to compute this to figure out these kokoro. Poke puzzles are sort of. They're like we look them up. Yeah sort of this cross between a su doku and a crossword. And you have these like sums of digits in the corners. And you're supposed to use that to figure out like fill in the digits in the middle of the puzzle. I've never played these but it's obviously like a nice little math puzzle that you can solve so he was writing a thing to solve these the Python version took PI point.
Five seconds to solve the original naive. Julia thing that he'd written was thirteen point. Four seconds and then with some minor tweaks it got down to like point nine seconds but then people kept sort of whittling away that was mostly just like don't you know the Python version was it was pushing stuff onto arrays where the other one was like creating a new ever a every loop. So you're like yeah okay allocating a new array. Every time is going to be slow but we can also push on to the ends of a race. So that's not a problem and then this this guy. Jason Merrill who's been around our community for a long time came like got into this and like really really optimized the heck out of it and he got it down to 13 microseconds which I was just like okay. Yeah we get like sometimes orders of magnitude performance improvements. But that's like ridiculous. That's six orders of magnitude. So how do you get a million. EXPEED up on this so the original version. I'm not going to show you the code for the original version. But you can imagine what it's doing it's like making a raise of digits and then sort of doing stuff to sum them. Whatever and the breakthrough that that the Jason had was he realized well. Hey the in this language it's actually really easy to make custom data types. That are extremely efficient but also convenient to work with so why don't. I just make a custom data type for this digit. Set thing because it's doing all these digits set operations so he creates a digit set type which is immutable which means it can't be changed. Which lets the compiler do nice things with it. It's also typical for like numerical things you want them to be immutable. Because they're sort of identified by their value not by some location in memory. There's sort of sort of like philosophy about that but I think the easiest way to see why immutability might be nice is that there was a bug in an early Fortran compiler where you could mutate the value of integers because it actually stored integers the same way it stored other variables so too was just like a name for a location in memory that happened to have the value to add it and the bug was that you could change it so you could write two equals three and then for the rest of the program to was three and like.
That's obviously crazy. You do not want to do that. You can have similar fun in. Python 2 by changing the values true and false that is allowed. They fix that in in. Python 3. So you know this is. There's there's certain things that really should be identified by their value if you change if you change the like imaginary part of a complex number you don't have that complex number didn't change. You have a different complex number right whereas with like a raise for example like you can have the notion of an array containing different stuff at different times is coherent that makes sense the array is like a container all right so we define this immutable digit set and the observation. Here is that to represent a set of digits 1 through 9. You only need 9 bits right and so an 16. Suffice a's that gives you plenty of room to work with different subsets of digits and then the rest of this stuff is just ways to construct a digit set from an array the iteration protocol so you know python has an iteration protocol you can make this thing iterable and Julia to the lot. The length so the length of a digit set is you count the number of ones in the bit representation of it and this is a one instruction thing on x86. I'll show you some some x86 code in a second. I know you're excited. And then how to show the thing so it's nice and usable and then some set operations so Union. Intersection sim diff and set diff and these are just defined in terms of ands and ORS of int of the int field. That's inside so let's evaluate this create one. I should actually. I should clear all the output to add suspense.
All right so there. We created a digit set. So you can see here that we can do things like oops the Union we can also use unicode union you can also write. Union but we the unicode union is an alias for the other thing so you can see the union of these two sets. You can see the intersection. No it's not. It's called cap inlay tech. These are the latex capes. So you see the intersection 356. I don't think that there is a fancy Unicode symbol for the symmetric difference. But that's the symmetric difference between the two. So here's what's nice. Check out the code for this ok so the native code for that is so that the first three instructions are actually just like stuff like set up for the function and the last two are just how you return from a function. The only thing that actually does anything here is this one or double or word instruction which is the x86 instruction for for oaring to 16 bits. Bit things together and you can see the same thing for the code. The code native for the intersection is just that same thing except with an and and the sim diff it's an X or so I mean this is pretty slick. This is very nice yeah so anyway also the count ones. Let's see what count ones i haven't actually tried this one but let's see so length of a is five ok code native. Oh no we blew it on that one. That's lousy codegen yeah can. I have to file a file an issue and fix that yes. Oh yeah okay well i mean we ought to be able to do better than this though ok so now now that we've got so the nice thing about this is you can see that you interact with this thing as though it was just this array but it's actually represented with this incredibly efficient encoding and all the operations are incredibly efficient so this is like this is one of the like the mantras of Julia's zero cost abstraction get a real nice abstraction but there is no cost to it whatsoever and that makes doing you know math stuff a lot a lot easier because you get the you want to think in one level and you want the computer operated a very different level and then you just need to connect those two so now you can just make this thing a lookup table because there really aren't that many of these combinations so here's the code to do that and then you can.
For example look up this guy and it is a 12 element array of digit sets. This is kind of cute. The size of s that array is 24 bytes so it's two bytes per and there's 12 of them so this is 24 bytes to represent each of these one of these things now doing the decomposition thing is just the lookup and so we benchmark it and so the first time you benchmark it the first time you run something it does jit compilation so remember. I said Julie is dynamic but it's compiled that means that the first time you run something you compile it so it's going to it's going to allocate some memory because the compiler runs it's going to take a little bit of time but you know you run things more than once the second time it's 10 10 microseconds to actually do this computation that was taking 13 seconds before. And it's you know there's some variation but that's it. Oh one thing. I wanted to show this Julia set code so the code native for this is just like a page full of instructions x86 instructions so that is this is pretty fast like and this is you know this is actually if you look at how this is defined so you know for example. Let's say we wanted to add two complex numbers so 1.5 plus 2m + 2.5 plus 3 m's so that's adding two complex numbers you know which of course works. Let's say you want to see which method does that so the the plus function is actually this generic function with 206 methods. We can look at what the methods are like. That takes a little while because there's a lot of them and you can see that this is all the ways of adding things that Julia knows about this would sort of be baked into the compiler in most systems but we actually just defined it in Julia. You can actually click on a link here and it takes you to the exact commit on github that you are currently running to the line where it was defined and you can see that the way it's defined is that you want to add two complex numbers you add the real parts and they're imaginary parts and you just make a new complex number out of it but that still boils down to that's like very high level of abstraction these are this is a parametric type you're like going through many layers of dispatch but it still boiled down to just like a couple of instructions this also this is what.
I was talking about with like breaking down barriers right like the social barrier like not only is it available but like you click a link and you can see exactly where its defined and you could change it if you wanted to. It would break stuff if you try changed it but like you know but you can. Yeah so anyway. Uh so. That's that's okay so now the final thing. I'm going to demo is this is going to get a little a little more technical but I kind of I think I think we can handle it. I think it'll be okay so this is an efficient implementation of Savitsky goal I smoothing. Has anyone heard of this. I had not heard of it until next last week but someone came up with. This is a really good example of this. We do a lot of stage programming in Julia. So macros are a form of stage. Programming jit is also a form of stage programming. We have this new thing called generated functions which is a like really unusual but really nice form of stage programming. Essentially lets you at the point. Where the compiler is figured out the types that it's going to call a function on but you can actually hook in and be like. Oh yeah for those types. I'm going to tell you what code to use for those particular combination of types and you can generate the code using Julia and then run it so I'm going to show you how this works so what we do is we. We define a type for this filter the Savitsky book the live filter. Oh I can't I get out of restart the colonel there.
We go alright so i define a type for this. It doesn't actually have any fields because the only useful information is the type parameters which are m and n m is the like the window the half size of the window oh that you're smoothing and n is how many moments you want to preserve so the thing about the civets giggle i filter is it. Preserves the moments of your data so the first central moment second central moment which is basically like the mean the variants and like then higher up would be like the skewness and the kurtosis so you can smooth the data while preserving those properties of it which is a pretty nice thing. This was like a paper when it came out an important paper that had never heard of until last week. So here's what we do. This is the actual function the code. That does this and this is going to be a little bit crazy looking. But it's actually doing this in any other system that i know of would be like 10 like a hundred times crazier so what this says is this says okay if you if if. I would like to generate a filter for this for a particular mnn on data of some type and smooth and smooth output. So we're going to write the we're going to smooth the data in data and put it into the array smooth which has been pre-allocated by someone else for us. And what you need to do to like to figure out how to do this. Smoothing is you actually need to construct a Jacobian matrix and then do a saul the backslash which is like linear least squares solution and then that tells you the coefficients for smoothing the coefficients for smoothing don't depend on the data though they just depend on like how many moments you want to preserve and how the wide your window is. So that's why this is a good target for this and essentially what we're going to do is we're going to generate like the optimal looped code that you could possibly write for this by hand for a particular m and n but we're going to do it as soon as you actually want to call that and I will show you that.
So this is the code that generates the code. That then runs which you'll see in a second. It's easier when you see the actual that runs. That's just some convenience methods the data here is just going to be a range from one. To a thousand it's not even like an actual materialized array. It's just a range object but it works. It doesn't matter we'll specialized for that so now we're going to do a three-wide filter which means that three before and three after is what we're going to use for our smoothing and the four means we want to preserve at least four central moments and we're going to apply it to the data and we'll see what happens here okay and it's taken a second so the reason it takes a second it's good as generating the code and then compiling that code if we do it again it's fast and you can see that it you know kind of kind of vaguely looks like the input data but it is you know it sort of flanges at the ends and it does some other stuff in the middle and actually so in this generation code. I did here right here. I save the expression that I generated into a global so that I could take a look at it later so what I did there is. I actually generated this code. So this is the code that was generated bat function. And this is what's actually running. And so you see we spliced in the actual coefficients that were computed based on m and N and so that this is the actually important part. This is the branch free middle of this that runs on most of the data. This is just a little fix-up at the beginning because you're actually running off the edge of the array and we just wrap around and the same thing at the bottom but you want to separate those out otherwise you you have branches in your inner loop and it's just going to be a nightmare for performance. This is the fast way to do it and so now we can see you can time this thing so the first time it takes a while. Because we're doing that code generation and then the second time it's 50 microseconds.
I haven't timed this against anything else. Because i don't i don't have on hand a good library. That does this sort of smoothing. But i think that's fast. I don't know sci-fi hasn't no okay. We'll have to do a comparison. See how that goes and the interesting thing here too is you can look at oh. I forgot to load stats base. There we go you can compare the moments and you can see that. The the first pair of moments you know so the left-hand side is the Act moments of the actual data the right hand side of the moments of the smooth data. You can see that you know negative 1 point 4 times 10 to the negative 15th is almost zero eight. You know 83,000 and then you know almost zero and then 1.25 times ten to the tenth so we are actually preserving pretty close the first four central moments and then we sort of like trail off and don't quite make it anymore which was not guaranteed so for was the what we actually wanted so the cool thing here is that like. So let's say. I want to do something with like you know a width of seven for smoothing and I want to do five central moments because like whatever I just call this and now you know it goes ahead and does that and we can see that the expression that it used was this other crazy expression here in most systems. I don't know how you would do this like I don't know how you would stage this programming like in c++. You could do this but you have to pre compile all the combinations that you might ever want which is a real pain because you might not know until runtime what like combinations you actually want or you know and what you can't do this is a library in C++ you have to do it as a template and then that's you know that means you need a C++ compiler to even do this in other systems. I think you could do this but it just won't be efficient so this I think gives you a nice marriage of efficiency and usability also the nice thing is that from the user perspective this is the library right respective.
I can see that this is like a little bit. You're like you could be looking at this and be like a little scary from a user's perspective. You don't need to know any of this like this thing just works. It's just a function and it's just fast and that's all so. I think that is that is a pretty. That's a pretty nice thing anyway. I think oh one thing I wanted to end with if you go to Julia box org all you need is a Google sign-in we could expand it to other things but that's what we got so far and you can just use Julia in a a for free in a notebook. I don't know why this is taking so long to load there. We go. And there's actually there's a there's a tutorial and you can just sort of evaluate stuff. It's notebook format so you don't really have to think much you just kind of let go through and evaluate stuff and then you can play with it so I have a bunch of stuff in there but you can start here with like you know. Start the tutorial and you can just run code so you know the first thing you do try is like one plus one equals two you know. R and you get a different random number every time you do this and then you can click through to like the basics so I encourage everybody to try this out. Give it a shot even if you're not interested in doing a project just yet give. Giuli a try and see how you like it. Thank you have time for questions. I me to take that as a yes. I mean I guess the closest would be like MATLAB or Ruby. It's definitely in the the algal school with like the end end blocks. Pascal also had similar syntax. It's nice when you're cutting and pasting things because that can get kind of gnarly in Python. Yeah it's when you so like it statically compiled languages. It's you compile everything and then you run and those are completely set cleanly separated in classic dynamic languages. There's no compilation so it's just all runtime. Stage programming is where you interleave compilation and execution in various ways. Yeah yeah yeah yeah you can.
There's a package called pie call. Which lets you maybe. I should demo it. I don't know uh III can't I'm not enough time to demo it but like yes you can do zero copy sharing of between numpy and Julia it the package was written by Steve Johnson who was also the author of fftw which probably people have heard of. It's like this insanely good. FFT library that you've all used Eve whether you know it or not. And he's just amazing and the stuff he did with. That package is just incredible. You can call back and forth for example between Julia and Python functions recursively and it just works like not recommended but it works so it's as good as see typically so to the extent that those systems are as good as see I guess they're comparable to the extent but they're not then it's not comparable. I think that you know there were very much trying to do. Similar things one of the nice things i think that. Julia gives you is that this is the standard semantics like you. Don't have to opt into a different compiler a different potentially different behavior like the normal behavior that everyone tested their program with is the one. That's fast of course. The downside is it's a different language so if you're really set on using Python those other things are going to do. That and julia is definitely not design. Answer your question okay. Yeah yes but it's not like Lua where it's like a real small component. That's going to be easy to embed it's like. I mean with this is like everything in the kitchen sink. We're working on that we're gonna we're stripping it down. And targeting embedded infra straight into in like embedded architectures. But for. Now it's like it's a ton of stuff but if you're willing to load a ton of stuff then yes you can't embed it cool. Great thank you much thank you.
Lots of outputs running simulations doing all sorts of stuff shuffle it around a cluster and that any of that stuff is a nightmare to do in any of these languages so I also had Ruby to tie it all together and I realized at this point that I had created a Rube Goldberg machine. If you're not familiar with Rube Goldberg. He was an artist. I think active in like the 50s and 60s and you draw these incredibly complicated machines. Where like you know this guy is eating and his spoon pulls a string which like flips a piece of toast up and then the parrot eats it and then that like makes the parrot go down and then it put like spills a bucket of water and yada yada yada and eventually like something useful happens. But you know it's just like cobbled together out of all of these crazy pieces and that's what i felt like i was doing now it wasn't using Python at the time had I been using Python or head Python been in the state that it is now back in two thousand nine maybe I would never start it on the Julia project. But but we'll see i'll talk a little bit about that in a bit so this is actually the for language problem but we sometimes talk about the to language problem and the to language problem it actually goes under another name that's older which is austere holtz dichotomy and this is after this was a term coined by John Ulster Holt. Who was actually the creator of tickle. Then he was familiar with that programming language and he observed at some point in like the 80s. I think that there seemed to be to like classes of languages. That were quite difficult are quite different and that were sort of. You know what things gravitated towards and they were systems languages and scripting languages systems languages are statically. Typed typically scripting languages. I Dana dynamically typed systems. Languages are compiled scripting. Languages are interpreted systems. Languages have user-defined types. Like you have a struct or something like that that you can define and then do stuff with scripting languages tend not to let you do that and instead they provide a couple of very rich standard data types that you can do all sorts of different things like link a dick right in Python.
It's this one like really really powerful data structure or pythons arrays which are actually the get called lists. I don't know why that seems like a mistake. That should've been corrected point but um you know you can use them as a cue or a DQ or stack or all sorts of like. It's a very powerful array data structure they also systems languages also tend to be fast whereas the scripting languages are slow systems. Languages are hard whereas scripting languages are easy so this is this is the dichotomy and this was sort of the state of affairs in like most most of like the 80s and 90s and now in the like to thousands. We're starting to see that dichotomy sort of come apart a little bit. We're seeing this this. The systems languages sort of take some of the the properties of the scripting languages and vice versa. So because of this dichotomy there's a standard arrangement that people have typically used to do things. Because you're like well you know. I want convenience. But i also want performance i want both of those things and the compromise is that for convenience. Use a scripting language at the high level. Something like python or our matlab but you do all of the hard stuff in the systems language so when you're using numpy for example like the the fast stuff the stuff that's actually doing stuff with arrays is all written in c or actually. Fortran right because you're calling gloss if you do a matrix multiply which is still written in Fortran. So that and that kind of works like this is a pretty good arrangement. This has been sort of the state of the art for 30 years or something like that. I have myself written many high-level wrappers to low-level functionality and then you know you're happy because you can you can get the performance but you get this nice ease of use to but it is not a perfect arrangement it's very pragmatic but it's got some issues so one issue is aren't the hard parts exactly where you could benefit from a better language right so like why are we doing the hard stuff in see that seemed like kind of perverse right like I actually did try to.
I was fixing us one of our benchmarks because I noticed that C was cheating and doing less work than all the other programming languages and took me like two hours to get this thing not to crash. I don't know what I was doing. I had forgotten how to do see index incorrectly and see and I was just like how god I'd forgotten how difficult this was. So so this is this is a real issue. This is you know. C programming is not super easy it also forces vectorization everywhere right. There's this sort of like you're trot pot. That like for loops are kind of immoral like it as it turns out there just like impractical for people it writing libraries and high-level languages. They're not actually a moral right. It's fine to write a for loop the only issue is it's going to be slow. If you're in one of these high-level languages you want to like push as much work as possible into the low-level language but sometimes it's really awkward and you end up writing some very strange convoluted code to try to do this. It also can often force like additional memory allocation right because you have to go through a lot of create a lot of temporary erase in this style whereas like in C or something you would just like use one array and then mutate it in place and that would be fine and you would not ever create any temporaries. There's another subtler issue that took me a long time to recognize which is that it creates a social barrier between users and developers right so in these two language systems the users use one language they use the high level language. The developers actually don't spend that much time programming in the eye level language.
They spend a lot of time programming the low-level language and use it being a user of the high level system doesn't make you better you know qualified in any way to work on the internals you kind of there's just a barrier there right and I even fight even though I'm like perfectly well according to my story not that good at C programming but like I've done a lot of c programming. I'm okay at it you know if I'm using Python or are I tend to like hit that sea wall and then I just stopped because it's too much of a pain to find out where the sea source code is and how things work and I just kind of even though I could. I just don't so maybe maybe that's the moral issue maybe I should but but people don't it's just not practical so Julie is a programming language tries to crush and shatter. This oyster holds dichotomy when. I we want to have our cake and eat it too. So you know these are the two columns and we do this. We sort of take mix-and-match pieces from both columns julia is dynamic. But it's compiled we let you have user-defined types and standard type. So the standard types are pretty powerful but you can also do your own user-defined types and I'll show you an example of that in a bit. It's fast and it's easy. How fast this is sort of the obligatory performance chart this is time on a bunch of sort of you know micro benchmarks. That aren't really super interesting. Except that they time things like how good are you. At recursion iteration scalar loops parsing integers printing stuff. You know doing a little bit of matrix multiplication. It's time relative to see the the y-axis is logarithmic. So you know being higher up is much worse and the languages across the bottom. Julia Fortran go JavaScript. Python Mathematica our matlab octave python is sort of where you like. That's as good as a good interpreted languages ought to be. That's sort of where you are but then you have these like compiled languages like Fortran and go and then julia is the odd man out because it's sort of its dynamic but it's actually within a factor of two of see consistently beaten people actually see the Sun real-life code to this isn't just contrived from Mike Reuben.
I think this is a much more interesting chart. This is the same y axis but the x-axis is actually how many lines of code normalized between whatever the minimum and the maximum one was for each benchmark each languages and you kind of see these two regimes. You have you have a long. You know you have this corner right here. These are the fast piled language. These are the system's languages. These are the scripting languages. So you have this kind of choice you like saw you can be up on this end of the curve or on this end of the curve. You can either be concise and slow or you can be fast and kind of verbose. Of course there's one point that's sort of in the corner here so that's where we want to be. I like I like being in that corner. So this is my data science stack today. Julia for numerical linear algebra statistical analysis. Yeah he's Julie for everything Python can give you a lot of this right like you could use Python for everything but the c line right so it gets you down to from four line which is 22 but to get all the way down to one language you need to do more and project sleep. I are sort of like pushing in that general direction so maybe in the future people will be using pi PI and like writing Python code and it'll be like fastest see so that's exciting. Um so for the rest of this talk. I am going to give you some examples. I'm going to do a lot of live. Coding and you know things might go horribly wrong hopefully they won't. I'm gonna demo how you can use julie as a general perfect language so not just a language for numerical computing. Although i'm going to give a numerical example but you know whatever you can also use it for data analysis the usual stuff like data frames and then plotting some stuff.
And then i'm going to show what hopefully is a compelling example of a superpower for a user and then Anna compelling example of how it gives superpowers to library writers too and so that's that's that I'm going to do in the Jupiter notebook so everybody familiar with the Jupiter notebook this is a Python audience so I assume so everybody knows about this right this used to be called the ipython notebook but now it supports like 40 languages or something like that as kernels. Julia I think was the first like real usable back a colonel aside from Python fernando pérez. Invited me and jeff and some of the other julia contribute contributors to come to berkeley and write this kernel and we went and did it was a fun week. Ok so this is this is what Julia code looks like you can do things you can say that max iterations is constant and then you define you can have a doc string it actually goes outside the function instead of inside that actually makes sense in. Julia but I won't get into why and you notice that I don't really put any types or anything on here. I just say you know the dude. This is not this is the Julia set that I'm going to compute. Doesn't really have anything to do with the name. But whatever it's a cute coincidence so you take two complex numbers. Z and see. And then what you do is you do this iteration step this Z you say Z equals Z squared plus C and. Keep doing that until it escapes. Beyond the radius of 2 and the number of tot like iterations. You have do before thing escapes. It gives you like a color which you then can visualize it and make pretty pictures which we're about to see so you just you know press shift enter and that evaluates the cell and you can see that you got a generic function now call. Julia with one method. And we can you know. Call it on one complex point and you can see. I'm a little too tall for this podium so I'm going to put this up here. We can see some other points like okay. That's a higher point. I don't know so takes more iterations to escape we can do this.
With a comprehension of two-dimensional array comprehension syntax. Here so this is similar to pythons one-dimensional array comprehension syntax but you can do arbitrarily many dimensions and get in this case a matrix so this is a matrix of the escape times across a certain range. I chose the range so that at least some of the numbers were interesting. And they weren't all twos. But you can you can see this is 501 by a thousand thousand and one alright. Well that's kind of a you know a lousy way to visualize things so let's actually load some packages that do more interesting visualization the first thing. I'm going to do is define a color map. Which is it's the color this. This package colors provides a couple of different built in color maps with this one's going to involve red and blue. I'm going to get a hundred points because max iterations was 100 so I'm going to go from 1 to 100 and then I'm going to use that color map to turn escape times. I see if you evaluate the color map. It prints it nicely. You can see the colors it's actually just an array of colors but the we know how to print that as a nice thing and you can see here now. I make an image out of it and you get this pretty picture. Can people see that reasonably well on the projector. I think okay. Cool one of the interesting things about the Julia set is that you change these numbers a little bit. You change the starting point a little bit and get kind of radically different things and you get a lot of different really interesting pretty pictures but you know going in there. Oh another interesting thing about it is that doing it that way is the Julia set if you take this. So this is essentially a four dimensional shape because it's too complex numbers as inputs which are each two dimensional. If we take the diagonal of this set it is the Mandelbrot set so that should probably be familiar to people so the Julia set and the Mandelbrot set are related in this kind of this kind of way but it's interesting to sort of explore this set so let's use this thing called interact which gives you.
You can take anything that so that expression. I I used to generate an image. I can actually just say put this macro in front of it at manipulate is a macro in Julia and give it a range for sliders and then we're actually going to be able to like you'll see in a second you can scrub the sliders around and get different images all right so now we've got that and we can change it. Which is this is a much nicer way to explore different things and you can see that it kind of gives you gives you different patterns and their / they're all quite beautiful so this is a little more sluggish than I would like so one of the things this is doing is this is generating a new array every single time and going through in you know and then constructing an image out of it this is sort of inefficient. It would be nice if I could just pre-allocate one chunk of memory and just redisplay that every single time and not have to like you know go through all this memory allocation because that's expensive so this is a version of that which is slightly more verbose but does exactly that so I use a let block to say that what these INR are so that I can actually use their length. That's just kind of nice so what I do is I create a data array and here you see one of the particularly Julian things. We don't shy away from talking about types so classically dynamic languages try not to talk about types. They sort of try to pretend they don't exist. We decided that you know the first thing. I mean the whole like the whole initial point of numb pile is to have like typed arrays in Python. Which is the Sun typed language and so in numerical computing inevitably. You end up having to talk about type so let's just have one good system for talking about types and make that universal so it's actually just baked into the language and the language knows about types so here we're constructing a two-dimensional array a matrix the element types are the elements are actually this.
RGB type with an unsigned fixed point. 8-bit value in there. But you don't have to know about but this is all sort of like custom types that are just defined in libraries but they actually have the exact representation you would like for the machine which is you can represent. RGB values is just like you know triples of aid but of 8 bits. Then you create this image and then you just sort of go over and instead of instead of constructing the thing every time we just take a for loop and go through and just update the image every time for loops are not evil in Julia. At least they're not going to be slow. They're actually going to be quite fast so now we can see that we can. We can scrub through this and it's significantly more responsive which is quite nice cool any questions about that so far all right cool all right. So fractals are cute but you know we're here for analyzing data not analyzing fractals so let's do a little bit of data analysis. I'm going to load this package called data frames which provides a data frame data type. You know very similar to the one that ships with our and very similar to the one that pandas provides so. I'm going to load this. Jared jared is here somewhere this is one of jared's data sets that i actually you know just played with at some point was like oh this is a nice data set. It's good for displaying them demonstrating things NYC housing it's got a bunch of columns about like different housing sales throughout the city. And you can see that the type of this thing is that at the data frame and then. I'm going to use this package. Called gadfly which is similar to the it. It's a grammar of graphics style visualization package in pure. Julia and when I say in pure Julia I mean like the whole thing is written in Julia there is no C code involved which is kind of nice. It means you can hack on it more easily but it from a user's perspective it kind of doesn't matter so ok.
The first obvious thing you want to do is. I'll get rid of the headers so that we can see a little more. I don't know plot some stuff so like square. Gross square foot versus estimated expense. That should have some sort of interesting relationship. Right when we look at this and it's kind of this like diagonal blob and you're like well yeah. I guess maybe that has a relationship so I don't know the first obvious thing to do here is like there are really big values and there are really small values and a lot of stuff concentrated in the little small very small value. So let's you know let's scale it to be log-log scale so this is this is how you add attributes. If you're familiar with ggplot this is sort of similar to the way ggplot works. Ok now. We're getting somewhere there you go it's you can. It's like pushed pulled everything out of that. Like far left corner. We'd we now actually have things kind of spread out in a nice clean way. Still not like a little bit would be nice two separate things out so one of the things you can do is you can color by one of the fields and I'm going to color by borough and you can see that now. Ok yeah color by borough this is getting a little interesting and we can zoom in a bit and you can see that yeah so. Manhattan is the blue up across the top and Manhattan is like significantly more expensive. You're sort of higher up for long this log log scale then you know say Brooklyn is next and then Queens and then the Bronx and Staten Island. I mean you're I bet a lot of people here from New York City this makes sense. I mean so you can do you can do other things. Pretty straightforward stuff like plotting a histogram of you know so since we we figured out that like you know its income and square footage thing or have this log log relationship you know if we take a relay to take a ratio that should actually give us give us some like fairly consistent. Stuff so we can take a histogram plotting by color of these things and you can see again like you know there's there's kind of two modes for Manhattan and then Brooklyn is sort of further down branches.
You know further down is this the same same data just sliced a different way. So yeah you can do data analysis pretty straightforwardly. I I think that's you know that's useful. It's obviously interesting and like necessary that you can do this but this is not one of the things that sets it apart because you can do this kind of data analysis in any of these tools so. I'm going to talk about some of the other things the things I think. Give you real like super powers so the next thing. I'm going to talk about this is super powers for users this example this guy. Patrick Usul dinger. I assume that as a German name he emailed because and we get this a lot this happens is pretty typical like I wrote Python code and I poured it to Julia and it wasn't faster. What am i doing wrong or you are a liar because you said this would be fast. Depends i prefer the you know the former but the other thing happens to and unfortunately our mailing lists are full of people who are like love hacking performance right so this is always like this is always kind of fun because everyone is just like. Oh yeah here are several like cool ways you can make this much faster and like frequently. It's like orders of magnitude faster by the time you get done so in this case it was an interesting journey and performance because they're a couple of like the original so the original. Python version took point five seconds to compute this to figure out these kokoro. Poke puzzles are sort of. They're like we look them up. Yeah sort of this cross between a su doku and a crossword. And you have these like sums of digits in the corners. And you're supposed to use that to figure out like fill in the digits in the middle of the puzzle. I've never played these but it's obviously like a nice little math puzzle that you can solve so he was writing a thing to solve these the Python version took PI point.
Five seconds to solve the original naive. Julia thing that he'd written was thirteen point. Four seconds and then with some minor tweaks it got down to like point nine seconds but then people kept sort of whittling away that was mostly just like don't you know the Python version was it was pushing stuff onto arrays where the other one was like creating a new ever a every loop. So you're like yeah okay allocating a new array. Every time is going to be slow but we can also push on to the ends of a race. So that's not a problem and then this this guy. Jason Merrill who's been around our community for a long time came like got into this and like really really optimized the heck out of it and he got it down to 13 microseconds which I was just like okay. Yeah we get like sometimes orders of magnitude performance improvements. But that's like ridiculous. That's six orders of magnitude. So how do you get a million. EXPEED up on this so the original version. I'm not going to show you the code for the original version. But you can imagine what it's doing it's like making a raise of digits and then sort of doing stuff to sum them. Whatever and the breakthrough that that the Jason had was he realized well. Hey the in this language it's actually really easy to make custom data types. That are extremely efficient but also convenient to work with so why don't. I just make a custom data type for this digit. Set thing because it's doing all these digits set operations so he creates a digit set type which is immutable which means it can't be changed. Which lets the compiler do nice things with it. It's also typical for like numerical things you want them to be immutable. Because they're sort of identified by their value not by some location in memory. There's sort of sort of like philosophy about that but I think the easiest way to see why immutability might be nice is that there was a bug in an early Fortran compiler where you could mutate the value of integers because it actually stored integers the same way it stored other variables so too was just like a name for a location in memory that happened to have the value to add it and the bug was that you could change it so you could write two equals three and then for the rest of the program to was three and like.
That's obviously crazy. You do not want to do that. You can have similar fun in. Python 2 by changing the values true and false that is allowed. They fix that in in. Python 3. So you know this is. There's there's certain things that really should be identified by their value if you change if you change the like imaginary part of a complex number you don't have that complex number didn't change. You have a different complex number right whereas with like a raise for example like you can have the notion of an array containing different stuff at different times is coherent that makes sense the array is like a container all right so we define this immutable digit set and the observation. Here is that to represent a set of digits 1 through 9. You only need 9 bits right and so an 16. Suffice a's that gives you plenty of room to work with different subsets of digits and then the rest of this stuff is just ways to construct a digit set from an array the iteration protocol so you know python has an iteration protocol you can make this thing iterable and Julia to the lot. The length so the length of a digit set is you count the number of ones in the bit representation of it and this is a one instruction thing on x86. I'll show you some some x86 code in a second. I know you're excited. And then how to show the thing so it's nice and usable and then some set operations so Union. Intersection sim diff and set diff and these are just defined in terms of ands and ORS of int of the int field. That's inside so let's evaluate this create one. I should actually. I should clear all the output to add suspense.
All right so there. We created a digit set. So you can see here that we can do things like oops the Union we can also use unicode union you can also write. Union but we the unicode union is an alias for the other thing so you can see the union of these two sets. You can see the intersection. No it's not. It's called cap inlay tech. These are the latex capes. So you see the intersection 356. I don't think that there is a fancy Unicode symbol for the symmetric difference. But that's the symmetric difference between the two. So here's what's nice. Check out the code for this ok so the native code for that is so that the first three instructions are actually just like stuff like set up for the function and the last two are just how you return from a function. The only thing that actually does anything here is this one or double or word instruction which is the x86 instruction for for oaring to 16 bits. Bit things together and you can see the same thing for the code. The code native for the intersection is just that same thing except with an and and the sim diff it's an X or so I mean this is pretty slick. This is very nice yeah so anyway also the count ones. Let's see what count ones i haven't actually tried this one but let's see so length of a is five ok code native. Oh no we blew it on that one. That's lousy codegen yeah can. I have to file a file an issue and fix that yes. Oh yeah okay well i mean we ought to be able to do better than this though ok so now now that we've got so the nice thing about this is you can see that you interact with this thing as though it was just this array but it's actually represented with this incredibly efficient encoding and all the operations are incredibly efficient so this is like this is one of the like the mantras of Julia's zero cost abstraction get a real nice abstraction but there is no cost to it whatsoever and that makes doing you know math stuff a lot a lot easier because you get the you want to think in one level and you want the computer operated a very different level and then you just need to connect those two so now you can just make this thing a lookup table because there really aren't that many of these combinations so here's the code to do that and then you can.
For example look up this guy and it is a 12 element array of digit sets. This is kind of cute. The size of s that array is 24 bytes so it's two bytes per and there's 12 of them so this is 24 bytes to represent each of these one of these things now doing the decomposition thing is just the lookup and so we benchmark it and so the first time you benchmark it the first time you run something it does jit compilation so remember. I said Julie is dynamic but it's compiled that means that the first time you run something you compile it so it's going to it's going to allocate some memory because the compiler runs it's going to take a little bit of time but you know you run things more than once the second time it's 10 10 microseconds to actually do this computation that was taking 13 seconds before. And it's you know there's some variation but that's it. Oh one thing. I wanted to show this Julia set code so the code native for this is just like a page full of instructions x86 instructions so that is this is pretty fast like and this is you know this is actually if you look at how this is defined so you know for example. Let's say we wanted to add two complex numbers so 1.5 plus 2m + 2.5 plus 3 m's so that's adding two complex numbers you know which of course works. Let's say you want to see which method does that so the the plus function is actually this generic function with 206 methods. We can look at what the methods are like. That takes a little while because there's a lot of them and you can see that this is all the ways of adding things that Julia knows about this would sort of be baked into the compiler in most systems but we actually just defined it in Julia. You can actually click on a link here and it takes you to the exact commit on github that you are currently running to the line where it was defined and you can see that the way it's defined is that you want to add two complex numbers you add the real parts and they're imaginary parts and you just make a new complex number out of it but that still boils down to that's like very high level of abstraction these are this is a parametric type you're like going through many layers of dispatch but it still boiled down to just like a couple of instructions this also this is what.
I was talking about with like breaking down barriers right like the social barrier like not only is it available but like you click a link and you can see exactly where its defined and you could change it if you wanted to. It would break stuff if you try changed it but like you know but you can. Yeah so anyway. Uh so. That's that's okay so now the final thing. I'm going to demo is this is going to get a little a little more technical but I kind of I think I think we can handle it. I think it'll be okay so this is an efficient implementation of Savitsky goal I smoothing. Has anyone heard of this. I had not heard of it until next last week but someone came up with. This is a really good example of this. We do a lot of stage programming in Julia. So macros are a form of stage. Programming jit is also a form of stage programming. We have this new thing called generated functions which is a like really unusual but really nice form of stage programming. Essentially lets you at the point. Where the compiler is figured out the types that it's going to call a function on but you can actually hook in and be like. Oh yeah for those types. I'm going to tell you what code to use for those particular combination of types and you can generate the code using Julia and then run it so I'm going to show you how this works so what we do is we. We define a type for this filter the Savitsky book the live filter. Oh I can't I get out of restart the colonel there.
We go alright so i define a type for this. It doesn't actually have any fields because the only useful information is the type parameters which are m and n m is the like the window the half size of the window oh that you're smoothing and n is how many moments you want to preserve so the thing about the civets giggle i filter is it. Preserves the moments of your data so the first central moment second central moment which is basically like the mean the variants and like then higher up would be like the skewness and the kurtosis so you can smooth the data while preserving those properties of it which is a pretty nice thing. This was like a paper when it came out an important paper that had never heard of until last week. So here's what we do. This is the actual function the code. That does this and this is going to be a little bit crazy looking. But it's actually doing this in any other system that i know of would be like 10 like a hundred times crazier so what this says is this says okay if you if if. I would like to generate a filter for this for a particular mnn on data of some type and smooth and smooth output. So we're going to write the we're going to smooth the data in data and put it into the array smooth which has been pre-allocated by someone else for us. And what you need to do to like to figure out how to do this. Smoothing is you actually need to construct a Jacobian matrix and then do a saul the backslash which is like linear least squares solution and then that tells you the coefficients for smoothing the coefficients for smoothing don't depend on the data though they just depend on like how many moments you want to preserve and how the wide your window is. So that's why this is a good target for this and essentially what we're going to do is we're going to generate like the optimal looped code that you could possibly write for this by hand for a particular m and n but we're going to do it as soon as you actually want to call that and I will show you that.
So this is the code that generates the code. That then runs which you'll see in a second. It's easier when you see the actual that runs. That's just some convenience methods the data here is just going to be a range from one. To a thousand it's not even like an actual materialized array. It's just a range object but it works. It doesn't matter we'll specialized for that so now we're going to do a three-wide filter which means that three before and three after is what we're going to use for our smoothing and the four means we want to preserve at least four central moments and we're going to apply it to the data and we'll see what happens here okay and it's taken a second so the reason it takes a second it's good as generating the code and then compiling that code if we do it again it's fast and you can see that it you know kind of kind of vaguely looks like the input data but it is you know it sort of flanges at the ends and it does some other stuff in the middle and actually so in this generation code. I did here right here. I save the expression that I generated into a global so that I could take a look at it later so what I did there is. I actually generated this code. So this is the code that was generated bat function. And this is what's actually running. And so you see we spliced in the actual coefficients that were computed based on m and N and so that this is the actually important part. This is the branch free middle of this that runs on most of the data. This is just a little fix-up at the beginning because you're actually running off the edge of the array and we just wrap around and the same thing at the bottom but you want to separate those out otherwise you you have branches in your inner loop and it's just going to be a nightmare for performance. This is the fast way to do it and so now we can see you can time this thing so the first time it takes a while. Because we're doing that code generation and then the second time it's 50 microseconds.
I haven't timed this against anything else. Because i don't i don't have on hand a good library. That does this sort of smoothing. But i think that's fast. I don't know sci-fi hasn't no okay. We'll have to do a comparison. See how that goes and the interesting thing here too is you can look at oh. I forgot to load stats base. There we go you can compare the moments and you can see that. The the first pair of moments you know so the left-hand side is the Act moments of the actual data the right hand side of the moments of the smooth data. You can see that you know negative 1 point 4 times 10 to the negative 15th is almost zero eight. You know 83,000 and then you know almost zero and then 1.25 times ten to the tenth so we are actually preserving pretty close the first four central moments and then we sort of like trail off and don't quite make it anymore which was not guaranteed so for was the what we actually wanted so the cool thing here is that like. So let's say. I want to do something with like you know a width of seven for smoothing and I want to do five central moments because like whatever I just call this and now you know it goes ahead and does that and we can see that the expression that it used was this other crazy expression here in most systems. I don't know how you would do this like I don't know how you would stage this programming like in c++. You could do this but you have to pre compile all the combinations that you might ever want which is a real pain because you might not know until runtime what like combinations you actually want or you know and what you can't do this is a library in C++ you have to do it as a template and then that's you know that means you need a C++ compiler to even do this in other systems. I think you could do this but it just won't be efficient so this I think gives you a nice marriage of efficiency and usability also the nice thing is that from the user perspective this is the library right respective.
I can see that this is like a little bit. You're like you could be looking at this and be like a little scary from a user's perspective. You don't need to know any of this like this thing just works. It's just a function and it's just fast and that's all so. I think that is that is a pretty. That's a pretty nice thing anyway. I think oh one thing I wanted to end with if you go to Julia box org all you need is a Google sign-in we could expand it to other things but that's what we got so far and you can just use Julia in a a for free in a notebook. I don't know why this is taking so long to load there. We go. And there's actually there's a there's a tutorial and you can just sort of evaluate stuff. It's notebook format so you don't really have to think much you just kind of let go through and evaluate stuff and then you can play with it so I have a bunch of stuff in there but you can start here with like you know. Start the tutorial and you can just run code so you know the first thing you do try is like one plus one equals two you know. R and you get a different random number every time you do this and then you can click through to like the basics so I encourage everybody to try this out. Give it a shot even if you're not interested in doing a project just yet give. Giuli a try and see how you like it. Thank you have time for questions. I me to take that as a yes. I mean I guess the closest would be like MATLAB or Ruby. It's definitely in the the algal school with like the end end blocks. Pascal also had similar syntax. It's nice when you're cutting and pasting things because that can get kind of gnarly in Python. Yeah it's when you so like it statically compiled languages. It's you compile everything and then you run and those are completely set cleanly separated in classic dynamic languages. There's no compilation so it's just all runtime. Stage programming is where you interleave compilation and execution in various ways. Yeah yeah yeah yeah you can.
There's a package called pie call. Which lets you maybe. I should demo it. I don't know uh III can't I'm not enough time to demo it but like yes you can do zero copy sharing of between numpy and Julia it the package was written by Steve Johnson who was also the author of fftw which probably people have heard of. It's like this insanely good. FFT library that you've all used Eve whether you know it or not. And he's just amazing and the stuff he did with. That package is just incredible. You can call back and forth for example between Julia and Python functions recursively and it just works like not recommended but it works so it's as good as see typically so to the extent that those systems are as good as see I guess they're comparable to the extent but they're not then it's not comparable. I think that you know there were very much trying to do. Similar things one of the nice things i think that. Julia gives you is that this is the standard semantics like you. Don't have to opt into a different compiler a different potentially different behavior like the normal behavior that everyone tested their program with is the one. That's fast of course. The downside is it's a different language so if you're really set on using Python those other things are going to do. That and julia is definitely not design. Answer your question okay. Yeah yes but it's not like Lua where it's like a real small component. That's going to be easy to embed it's like. I mean with this is like everything in the kitchen sink. We're working on that we're gonna we're stripping it down. And targeting embedded infra straight into in like embedded architectures. But for. Now it's like it's a ton of stuff but if you're willing to load a ton of stuff then yes you can't embed it cool. Great thank you much thank you.