Seminar with Peter Fackeldey "Awkward Array: manipulating nested, variable-sized data..."

So, I'm really really delighted to welcome everyone who's here in person and everyone who's connecting online, the Berkeley Institute for Data Science. My name is Kirsty Whiter. I'm the executive director here. Um, one of the projects that we house and support and celebrate um inside of BIDS is the scientific Python project and community. who ran a developer summit in Seattle at the beginning of May. And I was delighted to get to meet Peter uh when we were up there. Peter is a scientist turned software developer, so PhD in particle physics and now thinking about solving problems with open-source software. Um we had a lot of fun looking at some GitHub statistics. He gets a A minus in the automated statistics from some app. Um, but I can promise you because of the quality of his work that he's about to tell you about, this is absolutely an A plus. So, thank you everyone um for coming along. This talk is going to be a little bit longer than our traditional um seminars. So, we're going to run kind of probably all the way up to the top of the hour and Peter and and I are here for a little bit longer afterwards. So, you are very welcome to drop when and leave the room whenever you need to. But you are also very very welcome to stay um ask some questions and we also invite comments discussion uh beyond just questions. So that's enough from me. Thank you so much Peter and over to you. Thanks. I hope you can hear me as well. Um yeah so thank you so much for the kind introduction and giving me the opportunity to speak here. Um I want to talk about the orray project uh today and uh first I want to start uh with a slide about me. So I'm currently a postto um at Princeton and I'm working in the RS project. I will explain this in the next slide. Um my scientific background is in high energy physics. So I'm a member of the US collaboration which is the second largest collaboration at CERN since 9 years now. Um and in case you haven't heard about CMS, Atlas and CMS were the two big experiments that found the expos. Uh my PhD thesis was about pair production. Um with the CMS experiment um yeah but during that time I also developed quite a software background. Um so I'm programming in Python since roughly 10 years and during the time of my PhD up to the end of my PhD I basically went from a a power user to a contributor uh of the psychot domain stack um to being a development maintainer. Um yeah and I personally have a passion for number crunching software. Um so if you want to talk to me about J feel free. Yeah. And a little bit about us. Um so as I mentioned Psychab and Iris before so Psychab is our community project our domain stage on GitHub where we connect all the um projects that that we develop and then there's also IRS which is a software institute funed by the National Science Foundation and you can see on the left on the bottom right many institutes all over the US um are are part of this. this great I want to give a little bit of an introduction about high energy physics because this is basically where came from. Um and in high energy physics we are interested in understanding the fundamental building blocks of our universe. And this is a lot about uh particle interactions and what we usually do for example at CERN or other experiments is we shoot particles against each other and take pictures and that's how we probe these interactions and um the analysis of these pictures is basically in a very simplified way done in three steps. So first uh we identify from these outcomes of the collision what is interesting then we build combinatorics. So we basically group together particles and see where they originally come from and then we reconstruct physical properties like an invarant mass. And on the uh left you can see a picture um from the omega minus discovery and you see a lot of white lines. These are all the particles that came out of this collision. Um, but as you can see, only the green ones are of interest. So, we need a way to isolate the interesting ones. And then you can also see that some of these particles decay into others. Like for example, this gamma here, this photon decays into E plus E minus. Um, so we but we only see E plus E minus in the detector. So what we do is basically build a combinatorics of both and reconstruct the problem and then we do this with all these decay products until we can reconstruct back the particle that we're interested in the omega minus and then we are for example interested in the mass or some other physical property. Uh they actually did this manually uh in the 60s. So they were looking at these 2D images uh and they analyzed by hand more than 100,000 of these pictures and if I got this correct this picture was found here at Berkeley so to discover the mines nowadays um the experiments are big uh yeah are much bigger um so I have this accelerator roughly 100 meters underground and this ring um has four main detectors. It's CMS, LHCb, Atlas, Alice and for example at CMS we measure proton collisions and now we have like a 3D image of all the tracks you can see on the top right but the analysis of these collisions are still the same. So we basically from all these yellow lines we're interested in which ones are of interest which ones are electrons for example imons and then we have to build combinations of those to find the original basically decay and then we reconstruct those crypto and this is what CMS and atlas did in 2012 uh for the discovery of the expos and here this is reconstructed the bar m of electrons or muons and electrons and you can see in red there's this excess of the HP um that matches the data from a computational point of view physics has um several challenges. So we have to deal with a lot of data um like our detectors they measure 40 million collisions per second um and each collision is of roughly one megaby of data. So we have a data rate of 40 terabytes per second. Um every event so every collision is independent of the next one. Um which gives us some nice factorization opportunity. Um every event produces a var number of particles. So if you collide let's say two protons there can be one electron or two or even more coming out. Um so we have to use this jacket data and um the check measure collection of information like momentum energy position. So we want to represent our data basically as arrays of stress and yeah as I said before u one thing that we do quite heavily is combining the products into their parents and this is very similar to joints in in SD. I want to focus on the first part of my talk about this lower part because that was really what a ray solves uh uh does. So you can represent and manipulate jacket data very efficiently. So the open rate project um this screenshot from pick up with me um was mainly developed by the Iris project and primarily by people who worked at Princeton in the Pixie uh department and the physics department. Um originally it was started by Jim Pasti who changed positions recently uh and angles uh Jana Henry Andre and me who now basically development mainly the port of course we have contributors from outside and also from other institutes of the research project to give you um a direct preview of how our library looks like. So um on the top you can see how you can construct an array. Oh sorry um so this array has for example an outer length of three and inside there are records that we call them records which are basically um yeah named fields. Um and you can see that for example this y here has first a length of one and here it has two elements and then it has three elements then four and then five. So this is basically um adding a jacket dimension in in an array. Um what we can do now in our array is basically using number like expressions to work with these arrays. So we can say okay from this array array I want to select all the fields so the field y and then uh on the last dimension I slice the first element away and then I take the square root of this. If you do this in equivalent Python, you would have to write several nested for loops and track the temporaries and intermediates. And so this is implemented using uh kernels and slicing syntax. Um so we can basically bring down the runtime the memory footprint by quite a lot by just comparing to run directly with primitive Python code. So in case you're not familiar with how jagged arrays um can be represented um on the left you can see how users would write or would basically define an output array or jet array um and on the right you can see the memory layout. So um let's take this example. So we have again three um outer elements and each of them have different lengths. So the first one has three entries, the second one has one entry and the last one has two. And uh we basically store uh the following information in memory. So we store a flat array of order contents. Um and then we store also starts and stops for example. And the starts and stops we basically say where does the first uh three elements start and where do they stop? And with that information we can basically construct an output or generate. Now we can also do some manipulations for example as I show before removing uh the first element in the last dimension. So you can apply slicing syntax like with no paras. And um the output is then basically slicing weight to 50 slicing weight five slicing weight 40. And interestingly, you can do that by just adding a plus one to the start and just have a basically a shrink mu on the data. You can also do manipulation that is basically passed down just to the content. So it's not affecting the structure. For example, taking the square root of the energy. Um we still have three elements in the first event, one in the second and two in the third. So we don't have to manipulate starts and stops. The only thing that we can do is basically route the numpy split. But then there are also manipulations um that for example sum all the energies in in um that require for so that's why we have here the awkward do summation and the one stands for it's implied on the x as well. And what we do here is uh we're basically summing um the elements in the last dimension. So for example uh here we would sum 50 + 100 + 30 280 5 and 40 + 20 is 60 and um yeah what we do with the content is summing it up. So we have uh 180 50 and 60 as as the new content. Um but for memory layout we do not need starts and stops anymore because we reduce the jacket dimension and now it's a rectangular array. So we basically don't have to track starts and stops any so basically these three types of manipulations. So we have structural manipulations uh as I showed before this only needs starts and stops or revers offsets for example. Um and this is very efficient because we actually don't need data for this. We only know need to know the offsets for the starts and stops. So if we have some something that we need to load from disk, you can actually skip loading the data entirely because the slides don't need it. Um so this allows us to do some IO optimizations. Um data manipulations um is exactly vice versa. So it's independent of the structure. Um it only needs the flat contents of data. So what we can do here is we can route the past number nine kernels directly to the net data. But then there are also these more complex jacket manipulations where you need structural information but also need the data and for that we wrote dedicated kernels. Um there's a or a CPP library that is automatically shipped with our library. Then um as I showed you before we are often interested in combining elements from so basically calculating combinations of values inside our code array um and this is similar to join SQL so let's take a example here we have this array X and then we calculate the ex the combinations like all the two combinations so we would do one two one three and two three for the first one, the second is empty and then we have uh yeah quite some more for the last one and you can visually uh look at this basically calculating all the combinations here. Um by default like the combinations that are duplicated are skipped but they can be included of course and uh we can also do this across various axis. Then um we're also interested in doing a cartisian product between two arrays. For example, between all the muons and electrons in our um uh events and for this we have this uh function called cartisian because it's a cartigian product. Um and here we can basically do a cartigian product between x and y where you then do uh 1 a 2 a 3 a for the first element. The second is empty again. And then for the last one we have 4 C5 C 6 C 7. Then uh I also said in the beginning that we are usually interested in grouping um uh things semantically because uh this is just more comfortable for us physicists to work with. So what we do is basically we have this concept of record arrays where we group arrays into a meaningful structure. So we have this outer array of events and this one contains a jacket array of electrons and these may have physical properties like mass, energy, momentum and so on. So you would have to read this as something like events is a record array of particle types like neons or electrons. Uh and where each of these particles is a record different physical properties and these fields uh electron and mass um they are fixed length. They are typed and possibly named. So you can also have a record array that has no name. It's just basically a triple in high energy physics events typically look like this. So um one event is for example it has two electrons which each has transverse momentum pi which has positional coordinates um and then three neons with the same properties. Um but in the second event you may have uh one less sorry two less neurons but then you also have photons and these photons may have also different properties like it has some identification properties. Um so it's highly reg um it's not only we're not only having this raggedness in the number of basically particles but we also have different um fields. Yeah. And having these record arrays allows us uh to pretty generally u reconstruct physical properties because our array has this concept of behaviors and I showed you before that these fields can be named and if we know a name for example energy and momentum uh we can use this information oh this is energy and this is momentum to calculate the mass. So we can basically provide automatically properties um based on existing fields to calculate the multiply and uh for physics uh these quantities are abstracted into uh these behaviors via the vector operation. Um for more complex reconstructions a direct provides um uh two more things. So there's a lot of interoperability with JIT compiled for loop expressions with number, Julia and C++. Um and we do have various conversions to ML libraries. So you can convert to Python projects. I want to give you an example of how this works with number. Not sure if uh you're familiar with number but I hope this gives a good impression. So um we basically initialize an awkward array and what we're now interested in is uh summing over the field y in the last dimension. So we want to have basically an array that has the elements 10 0 25 21 because 1 + 1 + 2 + 1 + 2 + 3 is 10. So what we do uh in with array operations we basically say okay first we select the field y then we sum over the last axis and then we sum again over the last a basically eliminate the the x's and then we end up with an array of 10 0 25 and 21. um very familiar with array operations that might be intuitive but you can imagine this can become quite complex. Um so what you can do is also writing for loops just in Python. do we first allocate. So we have this function here. We first allocate an output array for zeros the outer length of the array. Then we loop over the records of this dictionary X and Y. And then inside of this loop we loop again over these elements of Y. Then we sum them up basically. And this is much more expressive and probably easier to read than the first base thinking of how the Lo works. Uh and if you do this you end up with the same output. The nice thing here is that you can decorate this with numberit and it gets just in time compiled and will be uh super fast. Um and this works because a registers it's its own types as no number types. Uh then we can also have a look at the performance and um if you do a time it on the time before this roughly takes 164 164 microconds. Um and when you run the D kernel it's it's much less more than a factor of 10. Um of course the second invocation here does not include the compiler. So for the first invocation of this kernel it takes a little bit little bit longer. Yeah. So, Aquar provides all these ingredients for energy physics that I mentioned before. Um, you can do slicing. So, you can isolate interesting tracks or particles. Um, it does support the combinatorics to find paricles. Um, and then it provides various ways to reconstruct physical properties. But it also provides uh much more. So, there are different array types, not just jagged arrays. You also heard about record arrays but there are some more. Um there are different back ends. Um we have quite a lot of high level map operations. I just showed you some but there are also like AK dome and so on. Um and there are other features for example there metadata handling XR named dimensions and so on. I first want to quickly go through different array types. So typically users interact with this highle array which is called ak.c capital array and I showed you examples before to initialize initialize this one and usually users would only interact with this array type because it is supposed to look at like a not which is something people are very familiar with. um it supports a pretty duplicator printing and inside of this printing it uh contains a type information that includes the shape and d type. So if you look at the right there's this show function and a type here it has an outer link of four then a variable length that's a jet dimension and then it has basically two fields x and y where x is of floor 64 and y has again a jacket dimension type in 64. Um yeah, similar to XR these high level arrays also do have this address something that XR um introduced to the arrays. Um yeah and what you can do is uh you can use no like indexing including fancy indexing other or cut arrays. So you can basically say uh I want to index this the original array this AK.num num is basically telling you how many elements are on the last axis and then I say exactly three and you can see here um that I'm basically slicing away all the y um fields that do not have um uh the number of elements equals to three. The next level down are layouts and these this is basically how we internally store our code arrays. So they have an XML like uh structure and um each of the different types of these layout arrays and uh they have different purposes. So for example there's a list offset array that one has offsets and contents and this the jagged list that I showed you before. Then there are record arrays and record arrays have fields. So there's a field X there's a field Y. Um, and then on the the very bottom of the spec, so to say, there are the lumpy arrays, which are just that bad data. And you can basically plug together any um of these types and construct arbitrary complex arrays. Um, and then the lowest level so-called forms. And this is very the concept is very similar to um it reminds me often of pi trees from from Jaw PyTorch for example um because here we can basically decompose um a layout into a blueprint which is just metadata and the actual arrays that are at the bottom of the leaves. And uh we can do this with two buffers. And we have this blueprint form um and the arrays. This blueprint form is basically just JSON like metadata to reconstruct the awkward array and um the arrays that we uh basically separated from the form is just a mapping from unique keys for example here node to data uh to an array. And if we reconstructed this chrome buffers, we're basically putting node to data at a point where you find the py node two and then we basically plug together again um given a form of the mapping to an awkward array. Yeah. So there are basically three levels of uh awkward arrays. That's high level the one that users interact with. Um there are layouts that implement all the jagged logic logic and instructs of arrays and stuff like this. Um and there are forms which basically separate the metadata structure from the actual um arrays and yeah the letter two are basically awkward internals. They are accessible but usually not needed to to work with. Um yeah the second part that array provides a different back ends. Um and so what's a backend for bucket array? It basically has two purposes. Um a backend manages memories for arrays and it also implements jacket kernels the ones that are not automatically there from from you can basically with one line squeet. So either in the constructor you can put back equal something or there's this function to change the back end of the existing array. At the default back end is uh called CPU. The memory of arrays is managed by numpy and kernels are basically taken from numpai as a CPG. Uh this one runs fully on CPU as I said it's default back end and awkward CDP implements additional kernels that basically do not exist because they have to work on checked uh arrays. The second back end um is cuda. So n is managed um by cup and um the cutters are implemented using pup existing kernels or we extended actually also with our own cooler kernels. This runs fully on GPU. Um yeah and as I said awkward provides check kernels um that are needed for for I showed you before an example with number for CPU case you can use number cuda to write custom that work together with this we also have a jacket uh this one can run on CPU or GPU um that's chosen by Jack because Jack manages the memory Um the kernels are implemented they use pulley J kernels and we have the possibility to do this because they are the so-called segment kernels. Um segment kernels are basically you have a reference array that's defining like segments and then you can say okay the first elements are all zeros and then it sums basically these all zeros together in in the in the original. Um we cannot just compile this because it actually requires us to know the shapes at compile time. This is something that cannot guarantee um but we can make use of operative. So we can have the gradients uh through a way using operative functions. Uh and then the last back end um this may be a bit confusing now but I hope I can clarify a bit later. So there's a thing called type tracer. Uh and basically we introduced this concept on our own. So there's this type tracer array that manages um the memory of arrays and there's uh the type tracer kernel library so to say. And the idea is that type tracer arrays they are pure poison objects. They do not hold any numerical values but instead they just um are basically duct typed arrays. So they look and behave like arrays but they don't have any values in there. Um and uh yeah the set module as I said implements basically a subset of all the array API methods that we need to work with array but yeah as I said more more on those later. Um the last thing I want to highlight here because that was also a big chunk of my work when I started at Prince is adding named dimensions and access and uh named access uh have several benefits. Um of course it's more readable uh and understandable and you can imagine that many supervisors in high energy physics reach for this so they can read their students code in a better way. uh they lead to less errors because you're not summing over wrong axis because you're consciously saying okay I want to sum over let's say all the electrons um and they can provide safety opportunities. So what we can do is if you do a binary operation that involve like additions um that involve few access with different names uh we can throw an error because it's likely not what was intended. Um, so for example, you should not probably want to add events and particles because it's just a nonsensic nonsense operation. Uh, there's also a nice write up by covered in the PT uh when named XR department. I can recommend to read this. It's basically motivating it from a computer vision point of view. Um, and yeah, doing a prototype in Python. Um but that was a big inspiration to add them to a ring. So the basics is basically you can add named access directly to a constructor of an or array. So you can so this array here has two dimensions and we can just use a tupil to annotate the first and second one. So the first would be then second particles. Uh you can also use a dictionary and basically map um a name to a dimension. So this would be the one means the first dimension, one means the second dimension. Uh or you can use this uh this function to add uh to an existing array dimensions. The nice thing about the dictionary notation is that you can also pin a name to like minus one. So you want to basically have the last dimension is always called articles like this. Um and then once they have these named X's attached you can use uh these names for for example for reduction. So you can use AK do sum on this array and you say I want to sum over the parts. Um so what we do here is basically one so 1 + 2 and 3 + 4 and we end up with an array three seven. And you can see in the representation it only has no events there because we reduce basically the particles dimension. It doesn't exist anymore. Um we do the fetch numpy operations directly. So you know instead of ak sum you can use np sum and will automatically translate correctly to the awkward one. And this one is also working um the name axis. Um if you sum for example over the events axis so basically you say 1 + 3 and 2 + 4 which is four comma 6 then uh the events axis is gone particles uh access is left and now particles points to zero because we reduce that direction. So it's not just uh the string is not just pointing to a fixed integer. it will basically dynamically adjust with the operations uh what the dimension means stands for um you can do uh of course you also have to adjust the indexes when you're doing indexing. So for example if you do select certain elements of an array oh sorry um you have to basically adjust the name dimensions as well. So for example, selecting um on the second axis the zero element means we do not have particles anymore. We just have events. If we select on the first dimension we do not have events but except we have particles again. Um if you add new dimensions like with non um you're basically introducing a new outer dimension in this case and this will lead to events and particle will be shifted by one because they dimension. So then it's going to belong to the second one as part. Um now that we have names we can also introduce a new syntax for indexing. Um so what we can do is uh we can use a dictionary for indexing. We can select with the name uh the dimension that um that you want to index and then here the zero or you can put a slice there. is the actual slice that is applied to the direction. Um, and you can also mix this with position and named axis. So for example, you can attach names to only certain axis and others are not having any names. So you can use a dictionary where using integer for the number of dimension and then the slice or the index that you want to apply. Um but you can also mix this or pause when you have um so here for example I'm taking from the first dimension the first element and then from the particles dimension I'm taking a slice of the basically zero to first element and yeah so if I do this here I'm basically selecting the one because I'm selecting the first group here the one comma rule and then I'm slicing zero to one which is basically retaining three length one dimension uh by just taking the element of one and as I said before uh we can also use some safety guards so let's cons these two arrays here um where I'm attaching um named ais of events and electrons and then the second array y I'm attaching events but now the second image has a different name um So if I add them um this will give a value error because the name axis are incompatible and it says it got electrons for the positional access of one uh which is a bit conservative as a policy but it's there for um basically providing some security and people can opt out of this. So once you attach an access you're not bound to use them. You can also for this alteration remove them and afterwards attach them again if you like. Yeah. So this basically um bring me to a little summary of array and this first big feature that I was working on. Um so I told you about the basics of oral arrays aboutics, the reconstruction of physical properties and various other features. Um like the levels of abstraction we have for arrays and different back ends. Uh and yeah in the last few slides I I showed you the new features. One thing I want to add for indexes is that um this is encoded throughout the whole code base of a brain. but also work with broadcasting um and various other things you can imagine in order access um and this is really a complete implementation of named dimensions or named axis. Now uh for the second part uh of my talk I want to talk a little bit about the first uh two or three bullet points that I skipped at the beginning that we have to deal with large amounts of data. Um and here we have to go to distributed computing and I want to highlight a bit how we facilitate this. So uh for high energy physics analysis we have um as a computational challenges we have a lot of data uh usually an analysis is working on roughly 10 to 100 terabytes uh I said before it's 40 terabyte per second of event rate or data rate but this is already for this 10 to 100 after the collaboration the collaborations do some central selection processing um And typically uh the collaboration provides files of roughly 2 GB of size compressed size. Um yeah so what we need to do is basically to scale to clusters that process these files or the data in parallel because uh yeah usually you cannot store 100 terabytes data on a laptop and even if you could then um you don't you do not have enough cost and the second part uh that I want to talk about in this scaling to distributed clusters is uh usually these files they are provided centrally from in the operations. They contain roughly one of 1,000 columns. So they have physical quantities of energy or mass of electrons. Um but typically physicists only need a very small subset of them like roughly order 100. Um and the memory of our clusters is quite limited. So would be quite bad if we load all the 1,00 columns but only use 100 because then we're wasting a factor of 10 um of memory. So we need a mechanism for loading only the required columns and this I call it IO optimization. This is something I want to talk about uh in a second this part. So now first about the scaling to a cluster. Um I want to start again thinking about uh how this works on a single worker on a single machine before scaling to whole field machines. So um we have I showed you before how an event in high energy physics looks like with electrons neons photons for example. Um what we're current what we're usually doing is so we group them in arrays like typically 100,000 because they roughly fit into a 2 GB regime and then that's what we call a chunk of events and this is usually the unit of data that's processed simultaneously on a single worker. Uh physics analysis looks like this roughly. So we basically apply vectorzed array operations from the original chunk or basically a process chunk. And here you can see uh event two is masked out because we're not interested in this one for example. And then the two boxes they annotate basically reconstructed some new properties that were not existing before. Um and then uh we're basically histogram uh for example this new blue um reconstructed uh uh information to a histogram and arrive at yeah for example x discovery and if you think about this this is basically a map reduce scheme so and there are solutions how to scale this to clusters and it's basic it's very simple so we can just process various chunks uh in parallel. So this map step is basically in parallel in a in a new dimension and then the reduction step is is happening again from all the different uh workers. And here we're basically scaling in this dimension with the number of workers and typically we're reducing the output on the client um in order to not have this 2 GB memory constraint of the cluster. uh we did implement a library called D awkward to make this uh automatically because here we want to connect D which is a P library to basically facilitate this separately uh and awkward array and D or allows you to build chunk or partitioned uh awkward collections. So basically um you work with a huge task or collection which internally has multiple partitions of chunks. Um and what it also allows you to do is build complex tasks. So not just map and reduce but also various granularity of uh the actual operations that happen. Now one example um sketch I want to show here when you go from the very top the original events to the histogram on the right. This is a small uh extra the actual sketch of the task. Um so these are just the two layers but they have more than few hundred ones. Um and here we have for example electrons or sorry letons on the very left. We select the first one. We do this combinations with the combinatorics that I show before. We have two electrons. we select them again. Then we access some properties of them and check if they are equal and stuff like this. So we can make uh each part of our analysis part of the hardware. Um but of course u wants to schedule each operation of the cluster. Um and we do not want to copy network all away all all the time the arrays between workers. So we basically got an optimization step which fuses together tasks to run on the on the same worker. Um but we also have the opportunity to eliminate tasks that are actually not needed. So if this is right um aate the mass of something but at the end we're not even involing it. So why in the first place? Um this optimization step can also get rid of unnecessary computations and this huge uh order of 100 layers graph is then typically optimized into the original use but it doesn't have any unnecessary computations for it. Um yeah so I think it's quite intuitive to think about this uh from two perspectives. So uh first from a single worker point of view um the single worker is just working on a chunk of events um we are using I call it implicit vectorization so basically sd for array paradigms uh multiers and so on to basically very fast implementations that work on these hundred thousands of events and the data locality is from the view of a single worker of course the data in the RAM have you um and and the CPU caches. From a cluster point of view, um basically a cluster works on multiple chunks of events at the same time uh each worker on one CH. This is the paralization that is also kind of implicit because we abstracted with D uh across a cluster of multiple workers. And the data locality here is that usually data is bringing via network into each worker and then it has it in RAM or CPU caches for processing. Uh yeah and now I want to talk about the second part just to remind you um here we have this problem of uh our files they have roughly 1,00 columns of different physical quantities but we're only interested in a subset of them. So we need to basically do some analysis of what we actually need and only load those into memory. And uh one thing I want to think you about is uh array tracing. And I hope you do not know the concept yet. So I would like to ask you what this function does. That's an idea. Adds one. It adds one. Well, uh this in this case it doesn't because uh you cannot write in Python like anything that if you add a one to it just does something completely different, right? Um so what we need to do is we need to find out what happens inside of this function and we need to understand what happens inside of this function uh first without loading any data and only what happens to arrays of this function and for this we use uh the concept of tracing um and how do tracers work. So this is a lightweight implementation of a tracer. So the tracer has this attribute which is a recording. Uh it looks and behaves like an array because it has shape and detail. And when you implement now the at functions we basically append to the recording the operation at and which arrays are involved. So sorry self and other and then we return a tracer with the new shape and d type. This is a simplifiedation. So in practice you would have to promote D type and so on. Uh and then if you feed this tracer with shape 10 and the type of flo we're basically adding to the global recording. Okay my computation needed uh the self and other. So in this case the array and the one um to run basically our array does not record which data uh sorry it does record which data touched or accessed unlike tracing injects for example does not record the computation itself because we're only interested in which data is needed to load it and not the actual operation. So that's why OKR implements its own tracer back end which is called type tracer. This is the back end I talked before and I hope this gets clearer where it's needed now. Um so the type tracer records only what data has been touched by the computation. It does not record the computation itself. And the tracing procedure is basically as follows. So at first we have to delay the actual computation. We don't want to run it because that means we can load everything. Um we can do this with Dask because Dask is when you write a code it just recalls everything in a computer graph and only uping comput uh it triggers the computation. Then instead of calling it uh with actual arrays the second step is we run the computation the analysis once but this time with a tracer not with actual data. Um and then we have this recording of all the columns that we need. So the first step is basically we not we load only these that are present in the in the recording from this and we do not load anything else. Um and then we run the actual computation but this time only with the loaded columns and ignore the rest and we made this automatically part of the D task graph optimization does awkward. So if you do code compute with awkward it will automatically try to do this optimal optimization for you. Um array tracing does have some drawbacks um because uh tracers do not know about their data. They only know about the shape and d type. So you cannot do any data dependent operations. So you cannot say if some value is larger than 10 in an array. Tracer doesn't understand it because it doesn't have any values. uh tracers can only follow one path of an if else statement and uh so in principle you have to trace both branches. Um and this is usually solved the jacks for example by using where statements. Um we need to trace uh the computation that would happen with actual data. So we need to make sure that uh we have pure functions. So if during tracing something you flip a switch and during the actual computation um it would do something different uh but we wouldn't have traced that path. So there would be a mismatch between what we trace and what's actually wrong. Uh and this is of course a bit tricky because it relies on the users to implement something like this. Um and of course it can be a little bit tricky to debug because um if you drop in like a debugger you sometimes you land in the tracing stuff sometimes in the actual computation and it's a bit confusing and uh turns out physicists are not necessarily familiar with this concept yet not as much as probably AI researchers. Uh so yeah it requires some support and teaching about these new concepts. So that's why we came up um with virtual arrays basically lazy arrays arrays as a second I optimization um I marked it currently in R&D it's available in the latest award release um and it's basically a concept that we've implemented a few years ago uh but we now revisit it um and the idea is basically we lazily load um the individual arrays of an awkward array whenever they are needed. So when like the sum needs some data or loads it um and what we did is we implemented this at the lowest model level form at the array container to have the highest granularity of basically laziness. Uh and the trick here is uh that each array in this container can be replaced by callable that would load the data. Um so every time operator needs it, it can just hold the recordable um and access the data. And this is you can imagine this is somewhat uh like generics map for for um and it turns out that liter arrays are much more robust than tracing. And this is because well if you do something like if this element is greater than 10 then well virtual race which is loaded at check um where tracing cannot because they never have access to loading data. Um it does have some drawbacks. So we're losing basically the ahead of time knowledge of all the columns that we need. So um this is giving us um not the uh ability for the most optimal data access. So we cannot move the request for multiple columns in a single network request and that just basically load them on demand and there can be uh some some performance hits if you will screen every time we are network just to give you an example. Just a time check because you're just at the top of the hour but you're really close to the end. Yes, I'm very close to the end. Um so just an example for virtual arrays. So this is basically um this from buffers that I showed before. So we have this blueprint form and uh arrays and uh this arrays container is having these unique keys mapping to arrays. But now you can basically substitute this by column list uh and basically for example putting lambdas before making even more complex like networks to to get the array. Um and then we get a virtual array. And what you can do now with virtual arrays um is you can look at them and you can see they have question marks here because they don't have loaded any data yet. They do know about their shape and D types because the form encodes all of this. If we access for example the Y field here, they still don't show us anything. But if we do a sum over the Y pairs, they load the data and call basically the callable from the method. So you have them available in memory. If you look again then on the right here we just loaded them. So we can have a look at them. And if you look at the original array you can see that well for y we have the values available but for x we do not because we have not triggered any computation of x. So there still is this question. So uh yeah overview about the distributed operate. Uh basically the upper level is called paralization. So on a single memory we use awkward array and SMD colors uh and then can also scale on the cluster with task and thus offers. Um for the IO optimization to reduce memory we can either do tracing to get all the required columns ahead of time um or we can do them on demand with good and that brings me also to my summary. So I hope I stayed within the hour. Um yes or is a generic array library for regular arrays. Um it focuses on working with regular arrays and not necessarily the serialization unlike p for example. Um it has quite some powerful features which are for example different back end metadata handling like xray and end access. Um andary can scale to large analysis needs. So we can parallelize on a cluster with task and we have this concept of IO optimization that lets you load only what you need and not. So thank you so much for your question. Thank you so much. Um it's amazing where so let me just come over. I have to stand near you to be near near the camera. Um so everyone online you're very welcome to uh stick around. I have a microphone so that the folks online can hear us uh when you ask any questions. So I want to just give everyone a moment to say uh if you if you do need to go if it's the top of the hour, if you've got other things that um are waiting for your time, thank you for coming. Uh but you are also very welcome to stay and we can ask some questions and um have a little bit of a discussion. So uh so folks of you are um online you can put a message in the chat or you can raise your hand and uh Lily is going to make sure that you um get the opportunity to to speak or if it's a very short question we'll just read it out. So does anyone have a question to get us started? Okay I see one hand raised. We've got a hand raised. Thank you, Matthew. Matthew. Yeah. Are you Do you want to unmute Matthew? Yeah, Matthew. Hello. You just voice it. Yeah. Uh, you know, great talk as always, Peter. Um so um uh I'm a fellow particle physics uh person so I'm very familiar with awkward but uh in your talk you focused a lot on the great technology that has been developed by the awkward team and also the applications to particle physics but since uh these like jagged and ragged structures are kind of generic can you uh maybe give some examples of like uh possible other areas of like domain science where this might be interesting or I know that you and Jim have given some tutorials where you've kind of given other uh examples as well. Uh yes so uh there are basically two application things that come to my mind directly. So the first one is text analysis or already actually implements uh uh functionalities for this. So um basically if you represented text uh sorry sentences as arrays where for example you have a jacket dimension because each sentence has a different number of words and then each word has a different number of characters. Um you immediately have a jacket jacket array. Then you can do things like I want to mask out certain words or slice the first word of each sentence and stuff like this. Um I think the second application that I know of is in uh astronomy. So you basically if you have um like uh yeah collection of objects in the universe like planets in a solar system um and you have the array of solar systems then you have multiple different number of arrays sorry planets per uh solar system. So this is also something where you can even represent a collection of product system to different number of planets. So um as I check the room thanks for the question Matthew any yes and just let me run to you with the microphone. Thank you. I can actually contribute another one because I remember I was on a project where you were interested in understanding for cars in a car company's car fleet and fleet as in they sell cars and all of those cars record data. think of Tesla or something like this um where you record data about each drive that that car would do um and record the battery discharge and then we were interested in modeling how the discharge behavior would affect the battery aging and that was also something that led us to I think we had like 15,000 cars and then each car had like different rides and that was like kind of um something that I remember we so that were XML data and it was like I don't know 10,000 XML files per day and then I remember we had to do a lot of processing on the cluster ending up with the same kind of raggedness. So um I imagine everything that involves this observational style of like real world objects that so yeah that's just like one an additional part here nothing special this sort of anything real world or anything real world kind of another question from online Daniel would you like to unmute Sure thanks um thank you I love the named axes feature that looks really nice. I've seen other takes on this and I I really like the take in awkward array. Are you part of the discussions to put named axes or dimensions in the array API? I think there might have been some chat about that in Seattle, but I I didn't catch it. And yes, so we had a discussion in Seattle about this. Um I think we uh did not follow up to you about this discussion. So, uh, it's it's the I would say open open to do. So, do do you mind just sort of um for the folks who went there just kind of expanding on? Yeah, sure. Uh so basically a few weeks ago in Seattle Python um community met um and talked about various things and um one of these were were uh specific for the array API which is still trying to generalize array libraries um in Python and uh I was uh bringing up the idea of adding name connections automatically to this. So yes, you can do it with master master race. Um I I talked with some people how did they did it for master race and I think there was constance of uh you probably can do this name access. Um yeah but unfortunately there wasn't a follow up about time during the workshop to to actually give it a try and implement. Thanks. And I I'll just I'll just add on top of that that kind of the one of the purposes of the scientific Python project and community is to recognize where um either there are multiple ways of doing something that makes it that are sort of that makes it difficult for users to kind of put stitch packages together or there's a sort of added complexity of thinking oh when I'm using this package I have to think about it this way when I'm using this package I have to do it a different way. So there's a piece of work around sort of harmonizing across the the scientific Python ecosystem. And then a sort of additional part on top of that is if someone's done something clever in one area, you can percolate that across to other other libraries and domains as well. So that's just painting that little extra bit of picture there. Do you have a question? Yeah, let me run to you. Wonder whether anyone used uh like for neural circuits because like you basically can simulate each neuron but the problem is that uh we probably don't have the the exact sort of a ground truth of each neuron. It's more like really just simulation because I'm just piggyback on his example use the flea the car fleet but um because I'm interesting working on the neuro circuit I'm just wondering whether someone has done something either directly using this or some somewhat similar I'm not aware of this um I would say I think the rule of thumb is basically if you have like some ragged dimension in your data like each neural has a different number of properties or uh values that you want to measure. Um then I would recommend to to use a regular array library. If you work with rectangular arrays anyway, um I think it's not necessarily um something you want to buy in unless you want to make use of the other features like different back ends, metadata handling and so on. Um I hope that you I'll also sort of editorialize on top as well which is did they oh maybe I can ask this as a question but I'll I'll throw why I'm asking it. um how young or old is awkward array as a project? And um part of the reason that we're really excited to have Peter here and to have you here is one of the challenges with um interdicciplinary research is actually having people from other domains know that this is a thing that you can do. So kind of being inspired by what's happening in high energy particle physics and translating that across to neuroscience or neuropformatics is a huge um goal and motivation of the existence of bids. And so I want to be very clear that just because people haven't done it yet doesn't mean that you should not be the first person to take it forward. So yeah, do you want to just answer the question of how how young or old is awkward array? Yes. So maybe Matthew can correct me because uh he's a bit longer than I accept than I am but uh from my uh knowledge Jim Kasi initially started a project roughly 10 years ago um and then it was first formulated in a zero as a first release and that was written in pure Python and numpy and then suddenly became uh obvious that there was a for these jagged kernels. So there was a rewrite for one to leave within C++. Um I think that was roughly 5 years ago. Um and then roughly two to three years ago there was another uh major rewrite. So there was a 2.0 which was basically still keeping um these kernels but much more focusing to uh the interoperable with the rest of the scientific Python ecosystem. So we can use task and checks as um yeah as libraries too. That's awesome. I'm gonna ask one more question but I know there's questions online. Sorry. Do you have an estimate of how many users there are for awkward array? I mean, you're welcome to just say a few some and lots, but you have uh I mean, okay, a little bit less than um number of users. That's a good question. Um I think Jim showed a graph about how many users there were over time. Nice. Um I do not have the number. Maybe he said many, many. Thank you. Thank you. You can say you had some good teamwork from me online online vende thousands of users that was expecting there was a question online chef on go for it. Hi Peter thanks for your talk really great to see what you've done with the library. Um, I was curious. So, now that you've got a labeled or named axes, um, it starts to look a little bit like awkward array can replace some X-ray functionality and you also get, um, acceleration free out of the box. So, so how do you see the relationship between awkward array and XR going forward? Yeah, that's good question. Um I think um okay I I might be wrong in this statement but I think Xray is not uh working with regular dimensions. Um it's basically an it lets you um basically partition rectangular arrays and then annotate um with metadata that is usually um focused on climate science um like let's say temperature um your arrays. So, awkward array does add this functionality for metadata handling. Um, but a much more focuses on the record arrays and working with record arrays where I had the feeling that X array is much more about representing partitioned arrays and working with them on the cluster. So, I rather say that last has much more similarity with XR than offer array. So, but I'm also not a big expert on XR. Yeah, it's Stefan because we we purposefully kept some time here for for discussion. So I do want to kind of invite like you know we don't we don't have to just spend an hour quizzing the so um I don't know Stefan if you have a thought about that because when I when I heard the question the way that Stefan um asked it I was like are you just like going to take over hex array does and actually honestly like that's a it's a slightly spicy question if that's the answer but it is also a feature to harmonize across the ecosystem if there's development that happens then it makes sense to bring things together. So, so is Stefan, do you want to unmute and kind of give your version of that answer or is it is it an open question for the community to friendly figure out? I was I was just thinking you know if I have like heavy rectangular array data I need to process maybe X-ray is still the place but it sounds like now if I have um you know an application where I would specifically benefit from acceleration for example and I don't need all of the X-ray features then maybe I can you know get get something out of um awkward array that I wouldn't have been able to do before. Anyway, just curiosity um I I think I imagine they serve distinct purposes. I was just curious um you know what what do you foresee users might do. Certainly I think I would agree with that. So there was one more comment just a comment that and I believe you had a question out the comment. Oh yeah, there was just a comment from Daniel. Um, he says, "As a happy user, both awkward and X-ray, I think X-ray does reflect some fine-tuning for the needs of the client science community, whereas Oh, then it just scrolled. Whereas awkward is a more fundamental data structure, right? I think one thing you may want to keep in mind is that uh um there might be cases when you use awkward array uh with rectangular arrays that you introduce uh starts and stops that are actually not necessary. Usually they're optimized away but maybe in some cases that might happen and then may blow a little bit of memory. It's not a lot but in that case it's say it's always safer if you have rectangular data. um cruise for example. Amazing. Thank you. And thank you Dan for um joining. Go for it. Yeah. So uh I'm I'm working on a project where we're currently trying to use awkward array for some optimization. Uh but we haven't been able to because I can't write from a Numba kernel to awkward array. Yes. Um and and I'm I'm barely started like looking into why that I think it more lies on the number side. But is there any like feature of awkward array that would be preventing it from having number kernels writing to awkward array? Uh yes. So awkward arrays are immutable and you cannot change values of an awkward array. Um but there's a solution for this with blockchain. So we have this concept of an array builder. um where you can basically write um or basically construct in for loops an or array um um with you do not use this AK array but instead of a different class um and you can populate basically your jacket array inside of a so that I think that may suit exactly your needs okay what would be the easiest these stories follow up with that question. Is that a question for all or is that a question for number one? You have an example on our documentation. I can see. Amazing. Amazing. Thank you. We are here to help problem solve. So, um I'm very interested in kind of the human side of sustaining total source projects. Um, what do you think is the most fun part of maintaining an open source project that's heavily used and what is the thing that keeps you awake at night worrying? Uh, I would say it's the same because uh because the fun part was about learning. You cannot just remove the feature because then a few hundred people will directly knock at your door and say uh what just happened? Yeah. Um this is uh fun to learn. I think this is a very fun experience. Um but at the same time uh this also a bit restricting because you cannot just uh overhaul the code base and introduce something new um that may break something. So it's always you have to think very carefully about uh what you change. Um so yeah, I would say it's actually yeah that's I love that. I love that answer. Thank you so much. Um folks online, thank you so much for for joining. Everyone here, you're very welcome to stay and informally um chat with Peter. Uh but otherwise, have a lovely afternoon and weekend when we get to it.