So, that gives us the SIR model. Now, we want to take Julia codes. We will do it in two steps. First, we'll just discuss Julia code for a single time step and then we'll put it in a loop that does all the time steps that we wanted it to do. So, this is the Julia code for a single time step. But now we want to discuss Julia functions in a bit more detail. We've seen a function briefly before, but now, let's discuss how a function gets the values that it uses for calculation. Here's an example and I will discuss it, but more fully. So, in this example, there's a line that just assigns the value four to b and then there's a line that creates a function. So, this line with this f and the parenthesis x, that on the left of an equation sign of an assignment means that we are creating a function which acts on x. What does the function do? It multiplies x by b and returns the value. So, if I've got b equals four and I've got x and now I want f of six. So, f of six means that x is replaced by six and we have b times six, but we know b is four, we expect to see 24. Indeed that is what we see. Now, why is that? If you are an experienced programmer, this will surprise you and that is because this value b is not actually inside the parenthesis. F is a function and x is the argument of the function. So, when the function is executed, we say f of six and so we know that the value of x must be six. But how does the function know what the value b is? The answer is that when a Julia function is compiled, it's compiled in a certain context and that is, any value that is visible in that context, and in particular, we are here inside this box. Actually all of the boxes are at the same value. So, all the values that are visible in all the previous boxes that we've executed will be visible to f of x when this function is compiled and so therefore we can change them. So, if we were to change four to a different value, let's say seven and recompile it and we run it again. So, now we've changed the value seven we recompile our function and we should get 42 which indeed is what we do. So, this is how we will pass the values of gamma and lambda to the updating function. So, when we write the function, we pass it only the values of S, I and R and the parameter values comes from the contexts. It's very convenient because it's parameter values are constant. So, every time we run the function, we use the same parameter values. So, you can declare them once and forget about them and that makes it really nice when we want to maybe experiment with different parameter values for the sake of fitting the model to the data because then we'll change them only in one place and we won't have to worry about parsing them to SIR in arguments. So, what we want to do is a function that updates S, I and R. So, we want a function that basically implements this model. It takes the current values of S, I and R and puts them in here on the right and it creates the new values S, I and R on the left. What we'll do is, we'll actually put S, I and R in a vector. So, the first element will be the susceptibles, the second element will be the infecteds and the third element will be the removeds. I could have called these S, I and R, but I wanted to just have the meaningful names. Now, we create the output values that we want, will be newS, newI and newR and they will go into this vector here. Note that, because we divide these values with spaces, they are rows and they will be just a single row. Eventually we might make many rows of a two-dimensional array, but this is just a two-dimensional array with one row and three columns. Note that this value is the output of this function. So, the return reserved word is very important one. If you leave it out, then the last of all the values computed by the function is what is returned. So, it's a good idea explicitly to put return in your functions you vary even slightly complicated because that gives you complete control of exactly what the function output is. So, now we just implement what we have. We have a lambda there, we have gam there. I use gam by the way instead of gamma because gamma is a built-in function and I don't want to overwrite it. Whatever the value of lambda and gam is in the context where this is compiled, that will go into the function and it will be used. Then the susceptibles is just the old susceptibles minus lambda, times the susceptibles, times the infecteds times dt. We implement our formulae and that gives us a function. One should always check a function a little bit. For professional code, one should do extensive checking. In fact, Julia is written where the tests of the functions are often written before the functions themselves and then the tests are run automatically. You should try and write your function in such a way that if it fails, it gives you some protection against the failure, it gives you maybe an error message or a warning. But we're not writing professional codes. So, we won't go into that, we'll just do a little bit of testing. So, we need a value for dt, we need a value for lambda, we need a value for gam, we need S, I and R values that we will put in the function and they must go in a vector. It's a nice way to do multiple assignment, you just separate your values on the right with commas and then for each value, you have a variable name on the left that matches it also separated by commas and that allows you to do multiple assignments on one line. Then we run the update SIR function on this little vector we've created and we see that if we start with 1,000, 10, 20, it updates to 975, 34.5 and 20.5. So, we see that the infecteds have gone up very fast, the removeds have gone up only very slightly and this they've gone up by 24.5. They've gone up by five and they've dropped by 25. So, the total is still 1,000. Okay. So that allows us to take a single step. So now we can go on to the loop structure. Now, the first thing is of course that the model can't be run indefinitely. There has to be a definite n time which we'll equal t sub f, and we will use t sub f equals 610, which is just slightly more than the number of epidemic days in the data from Wikipedia. So the number of steps that we need to take is in principle then tf over dt. But now let's just reflect exactly what is possible in terms of computers. It is not indefinitely precise to do computer arithmetic. So we cannot be quite sure that if we compute n like this, that it will satisfy this algebraic equation. We would assume, if we're doing pure algebra, that if this equation is true then this one is also true. But in computers, if this is our calculations not an equation, it's a calculation. Then we cannot actually be completely sure that this equation will also be satisfied. So the take-home message here is that in a computer you should not rely on exact equality even if you can reasonably expect it. What you can rely on is that a computer will be very close to the numbers that you expect, but they will not be exact. Well, it's possible of course to ensure exactness in certain professional contexts which require powerful computing and programming capacities, but we cannot rely on that. So we don't rely on exact equality. We only rely on closeness. If we know the value of n, if we compute it like this, then we can initialize an array to hold exactly all the values as they're collected. Of course if we take n steps, then there are n plus one values. So let's just write concept code to give us the loop structure. We'll use t final, as we said before, 610. That means that n steps will be t final over dt. But t final over dt is not necessarily an integer. So we have to round it off to be an integer. If you look up "help" on the round function, you'll see that you have to specify exactly what the type is that you want the output of a round to be. We just want it to be an int64, because that's the standard integer. Now we can initialize an undefined array. So the values in this array are undefined but they're all of type float64. So that's what these curly braces are doing here. It's an array with three columns and n steps plus one rows. Then we have the time values. Again, it's an array of float 64 values but now it's just a number of rows and we don't need to do anything more than specify that. So now that's the structuring which everything want to go, and in order to create it, we will specify initial values, we'll initialize the result values, and then we'll loop over exactly the number of steps and load the new SIR values into the time and into the results values at step plus one or rows. We'll also update the time and load it into the time vector. So this is in my opinion quite a good way to proceed. You go to the outer structure of what you want to do, you make sure that the logic is correct, that the loop structure is correct. It's easy to understand because there isn't a lot of detail, you just have to write a reminder to yourself of what should happen. I could have just made this one comment line you say, "update" and that's it. I could just have said, initialize over there. But in order to run the model, of course these details have to be put in and then we can run the model. So the full code now is all the parameter values we need so that lambda, that gamma, the time step length, we'll choose half the time step, that's the final, that's the s nought, that's the i nought, which must be nought zero, and there's the r nought, r nought isn't always zero of course. These are just thumb sack values that I'm just using to illustrate what's going on. They're not related to the Ebola virus disease at all. Then we initialize as we said over here, the initial values of s nought and of the result values. So we have to get the number of steps, then we can initialize the results, we can initialize the time vec. We first initialize results as, and completely undefined and then row one, all columns is s nought and i nought. Of course the initial value of the time vector is also zero. Then we can run it. So we update SIR on the vector that is in the current step value, from one end to the next one. So when step is one, this result value that goes in there and because step is one, one plus one is two, it will go into the second row. Similarly for the time vec is still step plus one plus dt. We go on to plot them. So again, it takes a while but you don't have to wait. So that's how the whole model code works, and now we can run it. So let's just run our model and it's very quick. It's already finished. We get the plots, I have already initialized plots in GR. So that's very quick. Now we can run the plot, it takes time. But just giving a plot the times, and the result values, we get three different curves as we want them. So that's very nice. The three columns, each column is plotted as a separate curve, and the legend is up in the right. Of course it doesn't read as well as it should. So let's just change the plot by adding a title and labels, and a legend. Now we can see that the susceptible start at 10,000 and they drop very quickly. The infected start at very near zero, they rise very very quickly. We see that here it's less than 50 days. So they peak, at about half of the total population are infected at the same time but then they die out very quickly. You can see here that the infected line is slightly below the susceptible line. So a few people were never infected. Similarly, the removes are slightly below the 10,000 initial value because a few people were never infected. The number of infected go down to zero of course because there's no more infections happening. So this is an illustration of how an SIR model produces output. Of course this is totally unrealistic for the Ebola virus disease numbers that we saw before. So what we need to do is to go on to try and adapt this model to give us the kind of output that was actually seen in Ebola, and then see whether when we discussed that, anything of insight, anything significant can be said about the Ebola virus disease epidemic. So that's what we're going to do in lecture two, and three, and four for the rest of this week. That's the end of lecture one. So that establishes the basics of an SIR model for the Ebola virus epidemic in West Africa. We are now ready to go on in the next lecture to do some interrogation of the model.