Teaching Numerical Methods with IPython Notebooks 2 | SciPy 2014 | Aron Ahmadia, David Ketcheson

David Ketcheson and Aron Ahmadia show how IPython Notebooks teach numerical methods, using a 2-D reaction-diffusion system as their example. They demonstrate interactive visualizations with Matplotlib, Plotly, Bokeh, and custom widgets that let students adjust parameters with sliders in real time. They discuss structuring notebook courses, handling version mismatches, and lab session logistics.
Chapters
Okay, so you should get this view. So put your sticky up if you are not at this notebook and you'd like help.
All right. Okay.
Taking over world plotting?
Here.
So my CSS file is actually lousy for the display up here, which is kind of interesting.
Try popping yourself in and out of full screen. I think that the browser's not redrawing the—
Right around. Yeah.
Okay, but you've all got it open on your machines as well.
Again, put up a sticky if you are not here, if you're not at this point. Okay, so one of the big advantages of the notebook, like I said at the beginning, is that you can include the code and the plots that are generated by the code, you know, side by side, and it can be edited by the student and they can play around with it, see what happens. You You can also put animations and interactive widgets. So that's what we're going to get into now. I'll go over very, very basic plotting stuff. This is not a Matplotlib tutorial. I think there was one yesterday. This will be extremely minimal where that's— with regards to that. But I will show some of the other options for getting sort of more than just a single static plot of your results. OK, so we're finally going to do a little bit of numerics here.
If you're not a mathematician or if you're not familiar with—
too familiar with partial differential equations, that's all right. But if you are, then you have some context for what we're going to look at. This is really just a backdrop for creating a plot.
So if it doesn't make any sense to you, just view it that way.
We're going to solve a 2-dimensional reaction-diffusion PDE, or a system of 2 PDEs as shown here. So we've got these 2 species that are diffusing. Importantly, they're diffusing at different rates, D1 and D2, and then they also react with each other. And this is a model that was proposed by Alan Turing back in the '50s for pattern formation, like how zebras get their stripes maybe, right? And you can actually generate lots of interesting patterns with these equations. We're gonna just use periodic boundary conditions. You We'll solve in a square box, make things simple. And then here are the particular reaction terms we're going to use.
All right. And here's the code that does this.
So you can see we're using NumPy. I don't think I'm actually using SciPy optimize.
That's probably a leftover import from some other version of the code.
I am going to use SciPy sparse because to apply a numerical discretization of the diffusion operator, I'm going to use a big sparse matrix.
OK. You don't need to worry too much about all this code because what I want to focus on is the plotting.
But if you want to understand where these plots came from later, you can come back and actually look through this code.
OK.
But basically we have functions that evaluate the reaction terms, a function that generates numerical discretization of the Laplacian, and then a function that takes a forward Euler step. So we're just going to use Euler's method.
We're going to use— yeah, we're just going to use Euler's method on this system of PDEs.
And we're discretizing in space with centered differences.
That's what this 5-point Laplacian sparse periodic does. It just gives us the 5-point stencil for the Laplacian. OK. So again, I don't want to spend time on that because we're already behind schedule. the focus here is plotting.
So as a first example, we're just going to create a plot of the solution of this PDE. We're going to run our simulation, create a plot using matplotlib. And the first thing you want to do is use this magic command, right? It begins with a percent, right? Aaron mentioned the commands that begin with percent are magic. They're called magic commands anyway. And what this does is instead of the plot appearing in a window, in a new window, In front of the notebook, it will actually appear in the notebook.
So that's what the inline option here is for. So do go ahead and run all of these cells or the ones below won't work.
So you can just kind of Shift+Enter through the thing. Okay, and as a first plot, we're just going to create this numerical discretization of the Laplacian, right? Our matrix A. And we're just going to use a grid with 4 points, and then I'm gonna spy that matrix. I told you I was going to use scipy.sparse because the reason for that is that this matrix is sparse. And you can see, so we have a 4 by 4 grid in space, so our matrix is 16 by 16, right? It has one row and one column for each unknown, and you can see that in each row there are only—
5 non-zero entries. Right, this is a sparse matrix.
But I did this just to show you, here's a matplotlib plot, right? We're not plotting the solution yet.
We're just looking at what the matrix looks like.
And hopefully you've used the spy command in MATLAB sometime before. It just shows you in a matrix which entries are not 0.
So all the white areas are zeros and the blue squares are not 0.
Okay, any questions up to here? Let's go ahead and solve the PDE and plot the solution. So I need to specify all those parameters that were in my PDE. And here I've pulled out of my hat some reasonable ones at the top there. And then I've got a function here that does all the sort of initialization. It creates a grid. We're just going to put random initial data. It turns out if you just seed this PDE with random initial data, it forms interesting patterns later on. And then I'm going to— here, these 3 lines will— this is matplotlib code to plot the solution. So I'm using this pcolor mesh command. If you've used pcolor in MATLAB or matplotlib, you know what this does. pcolor mesh is just a slightly faster version of pcolor with slightly less functionality. But for what we're doing, you know, plotting on a square grid, pcolor mesh will work fine. Yeah, this, this command will add a color bar. This image command says make sure the boundaries of the plot are tight to this square. I don't want to have white space hanging around. And also make the aspect ratio equal to 1, right?
Show it with the correct dimensions.
Okay.
Oh, and actually Aaron pointed out to me I don't need this anymore.
That last line of this cell is unnecessary, right? I think we'll find out.
Okay, so here's the actual code to run it. You can see there's a loop over time steps, and each time it calls this one-step function, which is just using Euler's method to advance our solution in time. I think I wrote at the top of the notebook, but I didn't say it, this is a crude method for solving this problem. It's not very accurate. It's not very efficient.
But for the purposes of today, it only requires a small amount of code. I didn't want to have a huge amount of code and make that part of it more confusing. Okay, so let's— okay, so I do need that show command for some reason.
Okay. You didn't run that command.
Did I not run the inline command?
The actual cell that runs it.
Oh, yes.
Okay.
Hahaha.
So notice all this cell does, all it contains are the parameter values here and then a function definition, right? And the rest of the cell is that function definition. So when I execute this, all it does is say, okay, now I know what that function is, but I didn't run the function.
All right. Thanks. All right.
And you'll notice something that we've added in the CSS, and I'll talk about it later, is that While the simulation is running, that cell is outlined in dark black. Normally in the default mode, the only indication you have that that cell hasn't finished running is the little asterisk in the brackets next to in.
All right.
Okay, and what's going on here? So I asked it to plot— I kind of did this a dumb way on purpose. So I asked it first to plot the initial condition up here, and then at the bottom there's some code to plot the final solution. Right, so I hoped to get 2 plots, right, and I only got one, except it has 2 color bars.
That's not very nice.
So what happened is that you— it only displays the last thing that was returned from your cell. Right, so it ignores the first plot. What I'm seeing is the final solution.
By the way, yours will look different because you had different random data at the beginning of it. It should look qualitatively similar, but Don't worry if the details aren't the same.
So it's just showing me the final solution. So there are ways to get multiple plots out of a single cell, like I could use subplots, I could explicitly tell it to output each of them, but I don't want to go into that because I'd rather do something more interesting. What I'd really like is to see the solution as it evolves in time because this is a time-dependent PDE. want to see how it goes from this random initial data to this pattern later on. And oh yeah, before I go on, another couple of things to notice that we'll try to remedy later. I can't zoom on this, right? If I try to drag it, I'm just dragging that image around.
It doesn't do anything.
I can't pan. I can't interact with it at all.
It's just a static image. I think it's just a PNG.
So I'd like to be able to zoom or pan in the plot. I'd like to be able to see its evolution in time. So in the rest of this notebook, we're going to look at ways to do those things. Okay, so this cell uses this IPython display module to repeatedly plot snapshots of the solution. So you'll see here it says if time is greater than next plot, right, so if it's time to produce a new plot, then We clear the output, we increase the time for the next plot, we sleep for a second just so that things actually appear, and then we display the new plot. All right, so if you run this one— have to run the next cell which actually contains it— okay, now you see the simulation actually evolving in time.
Right, okay.
And at the end I get 2 plots.
Now before I wanted 2 and I got 1, now I wanted 1 and I got 2.
You can do this more nicely, but it's a little more complicated to avoid having that second plot come out.
So, and I'm going to—
yeah.
Oh, you're right.
Oh, here?
Or I guess I maybe here, huh?
Close parentheses, plt.close parentheses all.
But I need—
yeah, okay, right.
plt.
Okay, thanks.
Okay, all right, and there's a nice—
there's a link here to more examples of using clear output with other types of plots.
Do you have a question, Lorena?
No, I was just excited. Oh, okay.
Okay, but there's a better way to do it, right, in my opinion. Okay, and So this is all just using built-in IPython stuff. Now I'm gonna go out a little bit on a limb and use some other packages that are all in beta.
Okay, so just keep that in mind.
Their interfaces are changing, their functionality is changing. So the first one is Plotly. So this is a web-based service. I run a command that actually sends a command to their server, they generate a plot, and they send it back, and it appears in my notebook.
Oh, that's cool.
Okay, it's free right now too. I don't think they have paid accounts yet. I think they plan to in the future make money, but I don't know. You know, I know, I know some of the developers are here so you can talk to them this week.
It's—
you can see that the code to generate this plot is very simple. And now this is taking a long time just because my internet connection is really, really slow.
Right, so that's obviously a caveat with this, with using this package. It's instantaneous with a good internet connection, and this is, yeah, really frustrating.
I can show you on nbviewer a notebook with Plotly plots.
Oh, except—
You have a bad internet connection.
Yeah, let's see, maybe it's the upload that's the problem. So this is actually running a complicated simulation. This is from my research.
Yeah, okay.
Here's a plot from Plotly. Here's the code that generated it. That's all of the code.
Right?
So not too much. And I actually created a layout thing here that includes 2 x and y axes and puts an inset in here.
Right?
So this is really fancy. But first let's just see if I click on this zoom button, then I can—
Zoom in.
Zoom in to part of the solution, right? See the details of it. If I double-click, then I get back out.
I also put a little inset here that shows the detail of the solution already. So that's another thing you can do with Plotly.
Is Plotly working?
Oh, there. Okay, it came up finally.
Probably we all hit Plotly at the same time, right?
And maybe crashed their server, and now they got it back up.
Actually, So there's something funny here, and I don't know what the right thing is, but you're all creating plots using my account when you run this notebook.
So now I have like 40 plots of these reaction-diffusion solutions.
You can click here and it will actually— it's creating a plot on their server, and you can go look at it on their server and you can modify it. Annoyingly, by default it shows up with a, you know, not correct aspect ratio.
In the notebook I got the aspect ratio correct because I set the width and height to be the same explicitly.
Oh, I see.
Okay, so you can do all kinds of plots with Plotly, and you can zoom and pan, right? Basic functionality that you would expect to have in a plot.
Can you zoom in?
Yes, so, oh yeah, and notice also as I mouse over, I see the simulation, the actual solution values. You asked if I could save the image. So if I click on it and go to the plot on Plotly, which I think I already opened over here.
Yeah.
Yeah, I can click save and edit here.
So I can probably right-click on it in the notebook and just say save image as, but I'm not sure. Oh wow, now we have a really serious aspect ratio.
Okay, so here's save, download.
Not sure what the difference— I guess—
oh, I can— yeah, so I can put a note on it, right? You can create a publication quality plot here interactively if you wanted, right?
Which is sometimes easier than doing it in a script, although less reproducible. But yeah, you can save it.
All right.
Yeah, so Plotly requires an internet connection. I know they're very interested in feedback on their product.
So if you use it, talk to them. They're here.
Another one, another package that's kind of moving in the same direction, but it's not a web service, it runs locally, is called Bokeh.
Hope I'm pronouncing that right.
I've— their interface is changing so quickly that I couldn't actually figure out how to get a pcolor plot working in the notebook. For this tutorial. So I just have a really, really stupidly simple sine wave plot, but it does show you what you can do with the Bokeh plot, right? So you can zoom the same way.
Okay, you can reset the view, you can pan, kind of pointless on this particular plot, but okay.
It has the basic controls that you would want to have, like you would have if you're matplotlib plot opened in a window, right, but it's in your browser.
Any questions about Bokeh or Plotly? I probably can't answer them anyway. Yeah?
It's probably just the internet connection. Is it still running?
Does it have an asterisk? No?
Oh, so you have to install it. Yeah, so open your— open a terminal and pip install Plotly.
Well, and then wait a while because your internet connection.
Okay.
It's a small package though. Okay, there are some other— I called them prospects. So other packages out there that are even more— maybe they're alpha, you know. So if you're really adventurous, you can try them out. There's mpld3, which is a combination of matplotlib and D3.js. D3.js I shouldn't try to define what it is, but it's a JavaScript library for interactive things on the web. So what this does is let you create matplotlib plots that then, you know, maybe— well, that you can interact with.
And there's a link to examples here.
Also Vincent and VisPy are 2 other packages that to me didn't look ready enough to show you here, but they aim to do the same thing. Or at least part of what they aim to do is interactive plots in the IPython notebook.
They may have other aims beyond that.
So just things to keep your eye on. And then what I've actually used quite a bit in the last year in my teaching is a little package called jsanimation. I think— so Jake VanderPlas wrote this. I don't think it's intended to be, you know, the final production version of this. My understanding is it's going to be superseded by something built on native IPython widgets. But it works very nicely with what it does. There's a little bit of discussion here of how to install it, which was on the Etherpad already. I think the easiest way is to just clone the jsanimation GitHub repo directly and then python setup.py install. But it doesn't— you can't pip install it, right?
So that's the level it's at, right? It's really like alpha. But it's works surprisingly well.
So here's just a setup that uses JSAnimation to create a playable animation of the whole simulation. And what I'm doing here, the only difference between this and what's above, is I'm creating this Python list called frames. I put the initial solution as the first entry in that list. And then every 5 time units, again, I append the current solution to to my list of frames.
Okay, so, so far nothing too exciting.
So it's running now and it's just saving all those frames. Okay, it's done. I didn't tell it to actually plot. And here's the code to actually plot it. It's a little bit complicated. So, so I'm importing IPython display from, from JSAnimation, right? So I'm telling it I want JSAnimation to do the displaying here. I create a figure. Right, and then I need to get a handle to a plot. So instead of using pcolorMesh, which would just generate the plot and return it, I'm using imshow here, which actually gives me back a handle to my plot. This was the trickiest part for me of using JSAnimation, right?
At least I couldn't get it to work using pcolorMesh. I have to use a command that returns a handle to a plot.
So this does that.
These arguments just tell it—
so here I'm doing an imshow on u, but I'm transposing it because the way imshow works it would show it transposed if I didn't.
And then I need to tell it the bounds of the plot, right, in both the color scale and x and y, and then how to interpolate between data points.
And then the most important part here, I have to have this function. I've called it fplot here. It's a function that plots one frame of my solution, and it takes as its argument the index of that frame. So all I do is I get my solution from my list of frames, and then I take my plot handle and I set the data to that new solution value, and I return that plot. And it should return the plot as a tuple, which is why this comma here is essential. And if this sounds annoying to you, Again, this is not a, you know, nicely packaged released code.
It's something that was put together very quickly.
Okay, and then here's the actual function call here. So this func_animation takes, well, 2— no, 3 arguments that are essential. First, I need a figure for it to plot in. Second, I need this fplot function that tells it how to plot one frame. And finally, a list of the arguments that will get sent to that fplot function.
One element of that list will get sent each time the fplot function gets called.
And then I can tell it how fast I want it to play, and there are a bunch of other optional arguments. But if I run this— so this is based on JavaScript, and here I get— okay, here's my initial data. That's the random data that we seeded it with. And I have these commands to play.
I'm gonna switch over to the play once bullet here, because I get annoyed by the looping, and hit play. Right, and you can see there, right, the evolution in time. Then I can go back to the beginning, play it again. Okay. Lots of fun.
So to me, this is the nicest of the options, right? With the clear output, we were able to play it once as it ran, but then we couldn't, you know, move back to a certain point in time.
With this, we can.
You can also do this with IPython widgets, which is what Brian is going to explain next.
Hi.
So I think that's it. And I'm way over time.
Any questions about JSAnimation?
OK. Can I say, test your microphone? Oh, look at that.
So we're really lucky to have Brian Granger here, who is one of the original authors both of IPython Notebook and the widget framework that he's going to show us in a sec. If you want to get to the link where he's at, I've updated the webpage. So if you do a refresh on our top-level webpage, you'll see that His notebook is first under Widgets Overview, and then the example that I had is second. So you want to use his example. Also, if you— he'll take you through the rest of it.
Yeah, I can show you how to find this.
You've got till 3:15.
Okay, great.
Thanks so much. If you want to go to nbviewer.ipython.org, I will show you where this notebook is as well as some of the other IPython tutorial materials that that we presented from yesterday. This also contains some other narrative documentation about the notebook. So to find this, the best way to find this is to go to nbviewer, and then nbviewer allows you to type a GitHub organization or user slash and then a repo name. And the organization for this is IPython and forward slash IPython In Depth. And so if we do that, It will take you to this page, and this is a repository that has all of these notebooks. The tutorial that we presented yesterday is Index Interactive Widgets, so I'll click on that. And then this is basically a longer 4 to 6 hour tutorial on both the IPython display architecture That was talked about previously, as well as the interactive widgets. And the notebook that I'm going to go through now is called Using Interact. So if you click on that one, and then if you want to, you can either follow along here in NBViewer, or if you want to download this and run it yourself, this button right here will download this notebook to your local system, and you can then run along. I should mention that this IPython In Depth repo also has in the examples directory documentation about the notebook, other aspects of the notebook, not just Interact and the widgets. So using Interact. Maybe take a step back for a second.
Sure.
So you've seen examples of situations where the IPython Notebook allows you to take content that's running in the Python kernel, bring it back to the front end or the web browser for visualization purposes. And we call that the IPython display architecture. And the main characteristic of the display architecture is that it's a one-way pipeline for data coming from the kernel, showing up in the front end, and then being displayed with a visual representation. And what distinguishes that from the widget architecture is that the widget architecture allows for 2-directional communication of data and visual representations from the front end back to the kernel, and then from the kernel to the front end. And that's based on a JSON message specification. Don't have time to go into that. What I'm going to be showing here is the highest layer in a multi-layered architecture that has sort of JSON and WebSockets and ZeroMQ on the lowest layer, an intermediate layer that we call the COM layer. Above that, there's actually widget objects, and then the highest layer is Interact, which I'll be showing here. Again, that longer tutorial that I showed you has more details about all of that. The reason I'm showing you this is that this API, the interact API, is where most users will want to begin when working with our interactive widgets. So here's some import statements. The module that contains these capabilities is ipython.html.widgets. The initial function we're going to look at is interact, so I'll run that. I should note that this does require IPython 2.0 or 2.1. This was not in IPython 1. So at the most basic level, Interact allows you to take a Python function and interact with it through user interface elements: sliders, checkboxes, dropdowns, text areas. And we build that user interface for exploring that function automatically. And let's look at how this would work for the simplest possible function that puts something on the screen. So f= print.
What's the output?
And if we pass that function as the first argument to interact, and then the second argument is a keyword argument, x equals 100, and the name of this keyword argument does need to match the name of the argument that we want to interact with in the function. So in this case, we're going to say x equals 100. We run that, and notice interact is able to infer from the fact that I passed it an integer, that I want an integer-valued slider for that user interface element. And when I drag this around, notice down here, every time there's a drag event, our architecture calls that function in the Python kernel and the result is sent back. If you pass a different data type to interact— so here in this example, we're passing a Boolean. Interact can infer what user interface element is appropriate for that data type. And so in this case, if it's a Boolean, the appropriate user interface element is that of a checkbox. Again, when I interact with that user interface element, the function is called on the backend. In this case, it just prints the value that's passed to the function. If you pass a string, the appropriate user interface element is a text area, right? An HTML text area. Notice when I call interact, I am following that command with a semicolon. And in Python, a semicolon hides the output of that command. Let's remove this and rerun it. And notice here the output that I was hiding is a textual representation of a function. What this means is that interact returns a function. In Python, what do we call something that is a function and returns a function? Anyone know this?
A generator?
It's a decorator. Python has this special syntax for functions that return functions. That special syntax is ampersand. Here is a function, ampersand interact, in this case is equivalent to g equals interact g and then the arguments. That's what that ampersand syntax is equivalent to. The reason that we return a function from interact is to allow interact to also be used as a decorator. Right, and so an alternate syntax that you can use for accomplishing the same thing is ampersand interact, pass the keyword arguments that you want to be used to construct the user interface elements, and then you've got the function.
Right.
And so we do this. So in this case, y is going to be a slider, x is a checkbox. That's why when you're using interact, you typically want to— when you're using interact as a function, not a decorator, you usually want to follow it with a semicolon to hide that sort of spurious textual representation of the function. Something else that comes up when you're using interact is that sometimes some of your arguments you don't want to have an automatically generated user interface element. A perfect example of this is if you wrote a function which takes something like a Pandas DataFrame or NumPy array as its first argument, And the other arguments of the function specify some transformation or some operation or some visualization to be performed on that single NumPy data array or Pandas DataFrame. In that case, you want that first argument that contains the data to be fixed. You don't want it to vary. And so we also have a fixed function that you can use to accomplish this. And I'll go up here to show you. So fixed is also imported from the same widgets module.
As interact.
Here's a new function, h , and in this case I'm saying I want p to be an integer-valued slider starting at 5, and q is going to be fixed at the value 20. When I run this and manipulate the slider, you can see down here that the second printed value here, which is q, is fixed. It doesn't change.
Right.
The idea here is that the keyword arguments passed to interact are a sort of widget abbreviation. They're a shorthand syntax that allows you to specify in a very concise manner what type of user interface element you want to generate for exploring that particular argument. We actually have a lower-level API that I won't have time to go into much here, but which we have a lot further documentation on is that we have actual Python widget objects. And so here's an example of one of these widgets, IntSliderWidget, right? So this is the type of Python object that's being built underneath the hood when you pass an integer value as a widget abbreviation to Interact. So you can specify the minimum value, the maximum value, the step, the value. And in fact, you can actually construct these widget instances by hand and pass them to Interact if you don't want to use the shorthands. And that's useful sometimes. Our shorthands have limited functionality. For example, if you want a slider and you want to specify maybe the description of the slider should be something other than just the name of the variable, you could do that. So here's an example of Interact where I'm manually passing an instantiated widget, and I get the same exact thing. The main thing that you have to know to use Interact is what are the shorthands for the different widget types. A few things about how this works. If the keyword argument to Interact is a widget instance with a value attribute, that widget is used as the user interface element for that function argument. Any widget with a value attribute can be used, even custom ones. So if you go out and you build a custom widget that's doing something else that our built-in widgets don't, as long as it has a value attribute, you can pass that into Interact and it will just work. Otherwise, the value is treated as a widget abbreviation that's converted to a widget using our internal mapping. That mapping looks as follows. So true or false is a checkbox widget. A string is a text area widget. A value, or a 2 or 3 tuple of integers, is an int slider widget. If those are floating point numbers, it's a floating point valued widget. A tuple of strings or a dictionary is a dropdown widget. And I'll actually show you examples of all of these. So here in this example, I'm passing a 2 tuple of integers. That's going to give me an integer-valued slider. So when I slide it around, the values in Python will always be integers. And the minimum value is 0, the maximum value is 4. If a third argument or a third element of that tuple is provided, in this case 2, that gives the step size, the granularity of the slider, in this case, It goes from 0 to 8 in steps of 2. The same is true if you provide floating point numbers. So here it goes from 0.0 to 10.0. When we run this, we have a floating point slider. Notice here the step size by default is 0.1. To modify that, again, we can provide a 3rd element of this tuple. So here we have 0.01. We run this. And we have the more granular slider that we wanted. You can also use these widget abbreviations in the decorator form of interact. So here I am saying go from 0.0 to 20.0, steps of 0.5. And notice in this case, I've already provided the widget abbreviation. Interact already knows how to construct the user interface element. The keyword argument to h itself does exactly what a keyword argument usually does, namely setting the default for that user interface element. Notice here it starts out— the slider starts out at a value of 5.5. It goes down to 0, which is this number here, goes up to 20, and the step size is 0.5. So with the shorthand notation and keyword arguments, you can get almost everything you would want to do in terms of customization of these slider widgets. Again, there are some customizations. For example, the description here, we're just pulling and inferring from the actual name of that keyword argument. The descriptions in IPython 2 cannot have LaTeX. If you do, I mean, it's fine, but it won't be rendered in the browser using MathJax as LaTeX. As of IPython 3, descriptions can have LaTeX in them and we will render that. Let's look at the case when the widget abbreviation is a tuple of strings. Notice we get a dropdown menu. The options available in the dropdown menu correspond exactly to the elements in that tuple that I passed. I can select oranges or I can select apples. Now, this particular shorthand always looks for a tuple of strings. You may think, well, what if I want a dropdown menu where the values are integers? You might think of saying, okay, that's nice. Well, then I'm just going to do 0, 1. What's the problem with that? It's not going to work. Am I going to get a dropdown?
No.
I'm going to get a slider, right? And so we have a different notation for getting a dropdown menu when you want the Python values to be something other than a string, and that is where the dictionary comes in. The idea here is that the keys in the dictionary will be used as the display names in the user interface. The values in the dictionary will be the Python objects that are passed to that function. When I run this, we see the keys here as the options in the user interface. When I select that, the thing that's printed, which is passed to the function, are the values of that dictionary. If you're using Python 3, Python 3 introduced the idea of function annotations. The syntax for that Is that immediately after the function argument you put a colon and then you can put any Python object. In this case, I'm putting true. This serves to provide inline the function or the widget abbreviation that will be used when that function is passed to Interact. I can't run this here. This particular example, I'm running Python 2. The main thing to note here is that this allows you to write library functions. So for example, SciPy could start to annotate its own functions with widget information that would allow any user to explore that function interactively by just passing it to interact without providing any additional information. If you're writing code that needs to work sort of both in Python 2 and 3, We do have an annotate decorator in IPython that can accomplish the same thing. It actually implements the Python 3 annotations in a Python 2 compatible manner. So we can run that, decorate that, and then when we run this we get the annotated version of the function without having to pass the actual widget abbreviation to interact. For simple functions like this, it's not too big of a deal to pass the keyword arguments to interact. However, you could imagine if it's a more complex function, something that takes 10 different arguments and the values, the particular way that you want to interact with that function is rather subtle, it could be hard to remember how exactly to pass all those things. You could bundle that in function annotation in this way. How are we doing on time? Want me to stop?
I think we can go into David's example.
OK. Great.
Thank you.
That was a guess.
All right.
Thanks, Brian. OK.
So just one real quick note since I imagine many of you ran into this. If you saw an error like this when you tried to run the Plotly example, that's because my code was written using Plotly 0.5.7, and they're now on 1.1.3 or something like that. But if you search their documentation on the web, it still is the old— the 0.5 API. So I couldn't find how to do it with the new one.
I found it.
You found it?
Okay, so We can—
I'll fix that later in the notebook so that it actually runs.
Apologies for that.
Okay, so I wanted to just show a couple of examples quickly of actually notebooks that I use for teaching. And the first one that I'm gonna go into, it's this aliasing.ipynb. Should be the first one in your list. Okay, so aliasing is this phenomenon when you try to represent a function, in the simplest case just a sine wave, on a grid, and that grid samples it at discrete points, right? And if you don't sample it at enough points, you can't really tell what frequency sine wave you have, right?
And if it's some more complicated function, right, then a more complicated version of that statement holds.
But essentially you need enough grid points to sort of resolve what the function is doing, and if you don't, then it will actually look like a different function, right? So So the word alias in English is like where something stands for something else. And in this case, you have something that looks like a low frequency but really should be a high frequency. So I actually used JSAnimation when I taught this in the spring in this notebook, but I rewrote it today as a Sage interact— or not Sage, IPython interact. Sorry, Sage had interacts even before IPython did. So I rewrote this as an IPython interact. Here's the code that does it. It's really simple, right? Most of it's just interact code. You can see here I'm plotting the function on a fine grid, and then I'm plotting it on this coarse grid where I've sampled it.
So let me just run this so that things appear.
Okay, so this first slider is the wavenumber of the function that I want to represent, and I've plotted the actual function in red. So if I change this, that red line changes. And then this second slider is the number of points at which I'm sampling my function. This is the number of grid points I have. So if I was solving, you know, some PDE, this would be the points where I'm actually keeping track of the solution. Okay, and I claim up here, right, that— or I claimed in what I was saying that sometimes if you don't have enough sample points, right, then a high frequency can actually look like a low frequency. Can somebody find an example of that and tell me the values of the wavenumber and grid points that—
19 and 10.
Sorry?
19 and 10. Okay, so let's go to 19 here. Let me just get it on the screen. It's a little bit hard to get these sliders exactly in the right place with the Touchpad. OK, yeah, so you can see— oh, does arrow— oh, that's nice. OK, thanks.
So you can move the sliders with the arrow key. All right, so we can see that this function is actually this very high frequency sine wave, but I can only see what one period— it looks like a function with just one period over the whole domain on this grid.
12 and 24. OK, let's go for 12 and 24. Okay. Yeah. Okay. So, okay.
So this is just an easy way for students to sort of understand this idea of aliasing. This actually becomes important— well, it's important in a lot of ways. First, this just puts a restriction on what kind of frequencies you can ever hope to accurately capture on a given grid when you're solving a PDE, for instance, or when you're doing other things where you have or you're sampling a function at certain points.
Also in Fourier pseudo-spectral methods, you can get generation of high frequencies that then get aliased down to lower frequencies on the grid, and this can lead to instabilities.
So this is an example of just a very short notebook showing one single concept, right?
And there's a couple of questions down here for the students to work on, right?
And this kind of thing, I might just have them look at this in class, right?
And then have them look at these questions and then maybe have just a quick quiz or ask a student to explain to the class what they've figured out. OK. Questions about this? Yeah.
What does it mean if you put the wavenumber to 13 and the grid points to 12?
OK. Ah, yes.
I was hoping somebody would find one of these.
Right? Yeah.
OK. So I mean, you can do the math and convince yourself. But it's sampling it just at all the points that are 0. So this is actually aliased to the 0 wavenumber, the lowest frequency mode.
Yeah.
So I'm not sure what you mean by what does it mean.
But yeah. OK. Yep. Very good. Yep.
This is also an example that's really good if you're teaching somebody statistics. or sampling theory, if they're sampling data from populations that have periodic data generation patterns, you can end up with data that looks like that when it's not.
Even though, yeah, okay.
Which we've used a couple times in stat classes.
Yeah, actually I was an astronomy undergrad and we studied variable stars, and you have to have a lot of observations with, you know, not too much spacing between them before you can be sure that you have a certain period and not some harmonic of that period. Or yeah, okay.
Yep.
Yeah, so here I'm creating a new plot each time rather than getting a line object and running setData.
Right.
I don't know which is better. In this case, they're both really fast. Any comments from the developers?
Here's the way to think about it. You have 2 different prices that you pay.
You have the price to pay to generate all of the data ahead of time, and then the only thing that you're paying after that is just the cost of grabbing it, right, from the code.
The second price that you So, given this thing, like, how long is it going to take? So, let me give you another example actually from the real world.
Real world. From broader general world.
If you're using Google Maps at maximum zoom, do you think all of those tiles have been generated ahead of time and rendered? Or do you think that they're being computed on the fly? Who thinks that they're being computed on the fly? Who thinks that they're rendered?
Rendered.
The most common ones are— they're all rendered ahead of time.
Okay.
Turns out that they're rendered and computed on the fly. All the really common ones on land are rendered, right? But if you're out in the middle of the ocean, out of nowhere, that image is actually computed on the fly for you. And so the right balance is somewhere in the midst of how much coding effort you have, how long it takes to compute that function, and how expensive it is to save those images.
So, but just to be clear about what the question was, you know, if you go back to the, the JS animation example at the end of this reaction-diffusion, right, here I was keeping my plot, keeping that handle to the plot, and just calling setData, you know, change the data that's in the plot that I already created. Whereas in this example, right, I'm creating a new plot every—
I'm creating a new plot on those axes every time. Right.
I could here get a handle to a line object and set data again.
Yeah. So okay, in this case I don't think it would ever matter because it's so fast. Thanks, Aaron. Okay, let's see.
We have a few more minutes. So there are 2 other much longer examples that are included in the notebooks in the repo for the tutorial. One is this multigrid and the other is this finite difference lab. Okay, MultiGrid is something that I would actually use during a lecture, but an interactive lecture where the students have their laptops with them and are kind of following along and maybe playing with the code. And then it has some exercises at the end for them to do typically as homework. This Finite Difference Lab is something that would be not so much a lecture as, you know, More a tutorial—
well, actually this isn't a good example—
but the kind of thing where students are primarily working. I'm not standing up at the front talking. I'm actually going around and just helping them answering questions all the time. So multigrid— so I'll go into this one. This is like sort of accompanying a lecture, right? Just showing the code and the results as you explain the concepts, right? So I would have the whiteboard next to me and be writing some of the math and then here running the code. Of course you could make it completely self-contained.
It just depends on your teaching style. You could have everything in the notebook.
So the idea here is to actually explain multigrid, but it's assumed here that the students have already read before coming to this the section in the textbook that I was using on multigrid. Okay, so how many of you know what multigrid is?
I don't know.
Okay, don't worry about it. It's a technique for solving systems of linear equations. And you use different— you use grids with different numbers of points.
So it's kind of related to the aliasing example that we just discussed.
I just realized I'm on the wrong—
Yeah, OK.
So if you— yeah, what you should do if you want to run this— sorry, this is pulled directly from my course materials. So in this— in the second code box that starts with %matplotlib inline, The last line of that, just delete clawpack.visclaw. So it should say, from jsanimation import IPython display.
There's a copy of jsanimation that lives in clawpack.visclaw, so if you have that installed, it'll work too.
Okay, so this uses jsanimation. I think I will, before I use it again, rewrite it purely in terms of widgets. But this is the version that I used in the spring when I didn't Maybe at the beginning of the course we didn't even have widgets.
I don't remember.
Okay.
So the idea is that it's explained a little piece at a time, right? So I give them a problem to solve, right? So we have this 2-point boundary value problem. Very simple, right? This is just 1-dimensional Poisson's equation. Here's our right-hand side. It's something just, you know, this function, all that matters is that it gives an interesting solution.
Other than that, you shouldn't care about it very much. And yeah, this is taken directly from the textbook actually.
So this is another way to use these notebooks as a sort of a supplement to a textbook that you're using. And here I'm just plotting the right-hand side function.
This is the forcing function in our boundary value problem. So I'm just using matplotlib there.
And then I'm going to talk about a way to solve this. And why is That's really strange. It must have something to do with the projector.
I don't know, but the bars in my fractions are disappearing.
Do they appear on your screens okay?
Yeah.
Yeah, okay. It's just something with the projector resolution here.
Okay, but I'm just using, you know, a centered difference approximation of the second derivative operator, right? And then I'm introducing Jacobi's method for this. So if you ever have taught a course in numerical methods, you're familiar with this, but it's an iterative method for solving the linear system of equations that arises through this centered finite difference discretization.
Right. And then in the—
so here the students actually are supposed to go implement Jacobi's method.
In this copy of it, I've just implemented it for you because we don't have time today. to go through the exercises.
So this just runs Jacobi's method 100 times. It just implements this equation right here and iterates 100 times.
Okay.
And then we're going to plot as a JS animation here the solution after every Jacobi iteration.
Okay.
In the interest of time—
We have a 30-minute break scheduled.
What?
We have a 30-minute break scheduled.
30-minute break scheduled?
Yeah, until 4.
Oh, okay.
So we still have some time.
Yeah, okay.
Are we okay going till 3:45 so we can get through this full example?
Yeah?
Okay, cool.
Okay, so what this is plotting, right, is the blue is the exact solution.
Right.
You can see our right-hand side f oscillates a lot more on the— our right-hand side of the equation oscillates a lot more on the right-hand side of the right-hand part of the domain.
And correspondingly, not surprisingly, you get more and more rapid oscillations over here in the solution.
This green line is our initial guess. So this is after 0 Jacobi iterations. And now if I hit play, it's going to show me what happened after 1, 2, 3, and so on Jacobi iterations up to 100. The red dots are the error.
So it's just the difference between the blue and the green. And what do you see?
So what part of the solution is the green line capturing well? The high frequencies. The high frequencies, which is kind of counterintuitive. You think the low frequency part is the easy part to get, but it's totally missing. The error now is very low frequency. So this is something that for me is really easy to see with this kind of animation. You see it just at the beginning there, very quickly grabs that high frequency information, but it just never irons out this big low-frequency error.
Right?
And so the idea of multigrid is really to somehow speed up the damping of that low-frequency error. And the reason that you don't get rid of that low-frequency error quickly with Jacobi is that you need— to get rid of low-frequency errors, you need communication of information between points that are far away from each other.
Right?
That's what low frequency is about. Okay, but Jacobi's method, if you look at this equation that we're solving each time, right, the update to u_i, hopefully— okay, so u_i means the value of the solution at grid point i, right? And then it only uses the values of the solution from the previous iterate at its neighbors to the left and right. It doesn't get anything from faraway values, right? So we have to take tons of iterations for the values in the middle to communicate with the values at the endpoints. And that's what we have to do in order to get rid of this low-frequency error. So what multigrid does is to actually use grids with fewer points in order to allow communication over large distances of that information.
That's one way to view it.
Another way to view it is that Jacobi—
you have here the average of the 2 neighboring points appearing on the right-hand side of Jacobi. Right?
So it's just sort of smoothing things out. And you can really think of it as doing that to the error. And so if you use a coarse grid, then when you're averaging over these large distances, you're getting rid of this low-frequency error.
Okay.
I'm giving a little bit of an example of how I would teach this here, right?
Okay.
But of course, really essential here is that the students would have coded up Jacobi themselves in class, right? Simple things like that.
I think help them get the concepts and internalize them a lot more.
Okay, so there's a discussion here about the convergence of Jacobi. So you can actually compute the eigenvalues and eigenvectors of the matrix that appears in Jacobi's method, and it turns out they have these eigenvalues. That is, the magnitude of the eigenvalues is magnitude of cosine of p pi h, where p goes from 0 up to—
2 pi h.
number of grid points you have minus 1, and h is the spacing between points. And essentially what this means— well, let me just show a plot here. So this is a plot of the eigenvalues of the matrix appearing in Jacobi's method. It's Jacobi's method when omega equals 1. So this is that absolute value of cosine. And what you see is that the eigenvalues in the middle, that is the middle frequencies, are getting damped a lot. This is what happens to each of those Fourier components of the error at each iteration. Whereas the components, the really low frequency components, their eigenvalues are close to 1.
So they're barely getting damped at each Jacobi iteration.
Also the really high frequencies are barely getting damped. But it happens that for the grid we chose and the function we chose, we don't have any frequencies out here in our problem.
Otherwise we'd see those sticking around too.
What it's really capturing well are these middle frequencies. But what we're going to do is introduce another parameter in Jacobi's method. So instead of just using this average, we're going to take a linear combination of that average and the previous iterate solution that we had. And the weighted average is going to be determined by the value of this omega parameter. You see I have 1 minus omega times my old solution plus omega times the new thing that I would have computed. So if omega is 1, then this just reduces to what I had before. If omega is 0, then nothing happens, right?
I just keep the solution that I had forever.
Okay, but this— if I change omega, then it changes the matrix that appears in Jacobi's method, and it changes how each of these things is damped. And the idea is that I'm going to use A value of omega right around 2/3 because this turns out if I want to damp just the upper half of the spectrum, all the higher frequencies, this is the best value you can use. And if you slide this around, you'll—
you can convince yourself of that.
If you want to minimize the absolute value of the eigenvalues in the right half of this plot, 2/3 is the best you can do.
Okay.
OK. Now the rest of this, you really need an hour and a half class session to explain.
So I'm not going to go into detail.
But the idea is that I'm going to apply Jacobi's method on a fine grid. And then I'm going to take away half of the grid points. So I now have only half as many and apply Jacobi's method on that grid. And then I'm going to do that on a grid with, again, half as many points until I get down to a grid with very, very few points and I can solve the problem exactly, right? And then I'm going to use what I did on all those coarse grids to remove that low-frequency error that I had in the initial problem. And yeah, I don't have enough time to explain it here, but I just want to go to the interact. So there are animations all along the way here showing the students the different parts of the of the algorithm. And if you just—
well, okay, you have to run all the boxes or the last ones won't run. So I think there's a discussion of aliasing here. So I pulled that out of here. All right. And then—
Okay, so there are some exercises here, and there's an interact that I coded up today. So again, I had done this all with jsAnimation. So this sets up a multigrid V-cycle, which is what's explained in the whole notebook above, and it lets you change a bunch of the parameters to see what happens. All right, so here, this is the solution that you get using multigrid, and this is the error.
You can see it's a lot better than what I was getting with just Jacobi's method, right? Now all that low-frequency error is gone.
So what these parameters mean, K is the number of coarser grids that I use, right? Or actually, K is the level of the finest grid that I use. Okay, so I'm using here 2 to the 8th minus 1 points because this is set to 8.
Okay.
Nu is the number of Jacobi iterations that I do to try to damp out the errors on each grid. And omega is this parameter that I explained before in under-relaxing Jacobi. And then there's also a Boolean variable here. I can choose whether or not to smooth after I interpolate. OK. And the idea is that for the students to understand this algorithm, they need to— if you just sort of throw it at them, they say, oh, well, it works.
But the idea is to understand how you came to that algorithm through a string of things that don't work very well.
And if you change these parameters— so if I don't use enough Jacobi iterations, say I just use 1 Jacobi iteration, on each grid, right, now the solution is terrible, right?
And you can—
and then I would talk with the students about why that is, right?
Have them figure out why do you need, you know, 3 Jacobi iterations to get a reasonably good solution.
Also, if we don't have a grid that, you know, resolves the solution well enough, right, then the solution will be very bad. And to some degree, interestingly, you can sort of compensate a little bit with, you using a finer grid and fewer Jacobi iterations or vice versa. You can even over-relax with omega here and see what happens there. So this is just more involved.
This would be like an hour and a half lecture where the students would here and there take 10 minutes to go code up a small part of it. But a lot of the code would be provided to them. And then exercises at the end. Okay, I think— well, okay.
Questions? I don't want to get too much into the pedagogy side, but if you have questions about it, I'm happy to give my opinion.
Everybody has a different style.
So, what are the prerequisites for this stuff? I found it difficult last year with a course that we did where we kind of The prerequisite was programming in something, and then we were using Python for the meat of the course. We didn't really spend enough time on the Python intro bit.
So what kind of course was this? Is this like a semester-length university course? And did the students program in class or only in the homework?
They programmed in class.
Okay. How many students did you have? Okay, and just one instructor?
3.
3?
Okay.
Okay, to me that sounds like a very doable— 3 instructors in parallel or sequential?
Yeah.
So, yeah. So I've taught this, this way with, you know, 10 to 20 students and a TA. And yeah, so let me back up. At the beginning of the semester, We have a 2-hour lab session that's just about how to use Python and the IPython notebook, which actually is still very short, right? And most of the students are coming from a background— these are graduate students, they have some programming experience, they programmed in MATLAB. I think for an undergraduate course, I don't know, it may be a little more difficult.
Okay.
So they've all programmed in MATLAB, but Usually none of them have programmed in Python, or very few. Yeah, so it's mostly learning Python syntax, you know, indentation, things like that, right?
The little things that trip you up.
And then most of the programming is really learned during the course, right? And it's through interaction with the instructor and with their peers. And that has worked fine for me. There may be a big gap between graduate and undergraduate.
I'm not sure. I haven't taught undergraduate students for 6 years.
So can you say specifically the problems that you encountered?
I don't know.
I think we didn't spend enough effort. You said you had a session at the beginning doing the Python programming. And they were doing different Python programs. I found that most of the students—
Yeah, so you go and you help them, right?
Yeah. Do you want to say more, Lorena?
No.
Okay.
Yeah, so the big advantage to doing this in class, right, is that there's someone there who knows and can help them instead of banging their head against the wall for you know, 3 hours and then learning to program it in the wrong way but something that works, right? You can actually help them get to the right way to program it.
You're grinning like you're—
Yeah.
Okay.
The what method?
I've never used that before. That's something Aaron has taught me, but I think I will in the future actually use Post-it notes in my class.
Yeah.
We'll talk about that in greater detail after the break.
I think that—
yeah, okay.
Students at the beginning of the semester, students have a really hard time asking for help, right, or telling you that their code isn't working or that they have no clue what to write, right? And you have to sort of force them to talk to you.
That's my experience.
Well, you just walk around. I mean, the way I use them is I don't lecture at all. And everybody's going through the material online in front of the screen. And I'm just walking around desk to desk. And the students, they don't have questions. But I have lots of questions. What are you doing? How does it look like? Why is it happening? Well, do you have any errors? And after a while, they start talking to you. And then they move on to the next guy. And you ask them questions. Do you understand that? What size is that array? Is that long or is it flat? You ask them questions.
So I should add that multigrid is probably the most complicated algorithm that I teach in this course. And so we've built up.
If you want to see a simpler, earlier in the course example, If you open this finite difference lab, you'll see that it's very, very basic and has them do some very simple things.
The other thing you do is have students interact with each other.
Absolutely.
Ask questions, surveys.
We've got one question over here.
Student response systems like Socrative. Ask questions and have them vote with their phones.
Yeah, so this could be the subject, of course, of another 4-hour tutorial.
Today we're kind of focused on the technological side. But yeah.
Well, I was just going to second that and say one of the good things about doing this is making the students code it right in front of you. Because I find that students will say they understand, but what they really mean is, I'm going to stay up all night reading the internet, and that's how I'm going to figure it out. She locked the door and said, no one's getting out of here. Wow. And that's how you really get the kids to want to come back.
Can somebody tweet that?
Lock the door.
So deep into the quarter, and it's really, really hard to get them up to speed. So yeah, you have to make it really fun.
That's awesome. All right.
Shut the door.
There you have it.
We, we are going to go to our break. We're going to come back at 4:00. Here's what's coming. We have 15 minutes on a little bit about how Software Carpentry teaches with the IPython Notebook and other general teaching techniques that we use. If you are a Software Carpentry instructor, feel free to commute with me and make sure that I'm not missing anything. From 4:15 to 4:30, David's going to talk a little bit more about the styling that we've been using, like David's style, his CSS style. And then from 4:30 to 5:00, actually From 4:30 till whenever, we're going to open it up. We're going to have a little bit more of a panel approach. We're going to take general questions and have general commentary. So if you have a lot of opinions, or you think that there's something that we should talk about, or anything along those lines, come talk to David and I during the break. And we can kind of make sure that that's sort of folded into the rest of the course. Thank you. We'll see you again at 4:00.
David. David.
More from David

HYDW08 | Prof. David Ketcheson | Solitary Water Waves Created by Variations in Bathymetry

Prof. David Ketcheson | Analysis and Modeling of Solitary Waves in Non-dispersive Models

Teaching Numerical Methods with IPython Notebooks 1 | SciPy 2014 | Aron Ahmadia, David Ketcheson

Teaching Numerical Methods with IPython Notebooks 3 | SciPy 2014 | Aron Ahmadia, David Ketcheson

Practical Experience in Teaching Numerical Methods with IPython Notebooks|SciPy 2014|David Ketcheson


