Hi Marco! Many thanks for this - great post! Just two quick questions - first, did you manage to find out where the problem lies in running

tcourse = ParallelTable[{sus[i, j], inf[i, j], rec[i, j]}, {i, 1, 500}, {j, 1, Mcities}]; // AbsoluteTiming

in Mathematica 10? Over here I am also experiencing the same problem, with it running unusually slow. Did you find a fix for that by any chance?

Second, I was wondering what would be the most efficient way to calculate the timecourse of the total number of infection cases over time in the last per-country simulation, and then optimize the model parameters against real longitudinal per-country data, e.g. as given in http://en.wikipedia.org/wiki/EbolavirusepidemicinWest_Africa (say for Liberia) Any thoughts perhaps? Maybe by running the simulation for a 3D grid of values for [Rho], [Lambda] and [Mu], and then using InterpolatingFunction to determine where it fits best? Or is there more efficient ways to do this?

Hi Vitaliy, Yes I had been trying that code too, although there too I had some issues towards the end of the code being very slow. In the code above I found out in the meantime that using

cities = CityData[{All, "Nigeria"}][[All, 2]];
names = cities[[All, 1]];
citypop = QuantityMagnitude[CityData[#, "Population"]] & /@ cities;
citycoords = CityData[#, "Coordinates"] & /@ cities;

or setting the global option

SetSystemOptions["DataOptions" -> "ReturnQuantities" -> False]

Solves the speed issue in Mathematica 10 - apparently working with Entities makes everything really slow. Using

EntityValue[#, "Name"] & /@ CityData[{All, "Nigeria"}]

strangely enough also turns out to be really slow. Well, any advice appreciated!

I'm completely new to Wolfram Language and had to sign up so I could comment and see if someone could solve my problem. I'm doing a maths expriment in which I graph the actual spread of Ebola and compare it to the theoretical SIR model. I found this and it's fantastic, but all I need is the SIR graph in the middle of the post and I'm not really sure how to get that. I'm also not quite sure how to manipulate Lambda and Rho values though... Can anyone help?

Hi,

In this post I was using the discrete version of the SIR model. In the video you linked they are using the continuous, i.e. ODE version. This could be implemented like so:

(*Clear everything*) ClearAll["Global`*"];

(*Define the equations*)
eqns = {
D[Sus[t], t] == -\[Rho] Sus[t] Inf[t], 
D[Inf[t], t] == +\[Rho] Sus[t] Inf[t] - \[Lambda] Inf[t], 
D[Rec[t], t] == \[Lambda] Inf[t], 
Sus[0] == 0.99, Inf[0] == 0.01, Rec[0] == 0};

(*Define the parameters*)
params = {\[Rho] -> 0.2, \[Lambda] -> 0.1};

(*Integrate the equations*)
sols = NDSolve[eqns /. params, {Sus[t], Inf[t], Rec[t]}, {t, 0, 100}];

(*Finally plot everything*)
Plot[{Sus[t] /. sols, Inf[t] /. sols, Rec[t] /. sols}, {t, 0, 100}]

One of the things one can try us to use ParametricNDSolveValue.

(*Clear everything*) ClearAll["Global`*"];

(*Define the equations*)
eqns = {
   D[Sus[t], t] == -\[Rho] Sus[t] Inf[t], 
   D[Inf[t], t] == +\[Rho] Sus[t] Inf[t] - \[Lambda] Inf[t], 
   D[Rec[t], t] == \[Lambda] Inf[t], 
   Sus[0] == 0.99, Inf[0] == 0.01, Rec[0] == 0};

(*Integrate the equations; rho and lambda remain as parameters.*)

sols = ParametricNDSolveValue[eqns, {Sus[t], Inf[t], Rec[t]}, {t, 0, 100}, { \[Rho], \[Lambda]}];


(*Finally plot everything*)
Manipulate[Plot[sols[\[Rho], \[Lambda]], {t, 0, 100}], {{\[Rho], 0.2}, 0, 0.3}, {{\[Lambda], 0.1}, 0, 0.3}]

Cheers,

Marco

That's fantastic! It worked, and well, so thank you!

Indeed, thanks Marco for your interesting work here, and also of course Vitaliy for pulling together that impressive post!

hI! In name of memory usage, the first part of the code I would suggest the use of Open read, as follows:

  SetDirectory["/home/etc...."];
    Airports =(*Import["/home/etc.../airports.dat","CSV"];*)
      OpenRead["/home/etc.../airports.dat", 
       BinaryFormat -> True];
    Routes =(*Import["home/etc.../routes.dat","CSV"];*)
      OpenRead["home/etc.../routes.dat",  BinaryFormat -> True];
    airports = 
      Flatten[StringSplit[#, ","] & /@ 
         StringSplit[ReadList[Airports, Record], "\n"], 
        1] /. {"\\N" -> "Non", "\"\"" -> "Non"};
    routes = StringSplit[ReadList[Routes, String], ",", Infinity];
    DumpSave["airports.mx", airports];(* here we save files in the format .mx*)
    DumpSave["routes.mx", routes];
    Close[Routes]; Close[Airports];
    Quit[];(*finish kernel*)

And then, for Airports:

SetDirectory["home/etc.../routes"];
<< airports.mx;(*!*)
<< routes.mx;(*!*)
GeoRegionValuePlot[
  Table[GeoPosition[ToExpression[airports[[i, {-5, -4}]]]] -> 1., {i, 1, 
    Length[airports]}], PlotStyle -> PointSize[0.003], PlotRange -> 1, 
  ImageSize -> Medium];
Quit[];

Routes : I Use "linksclean = Developer`ToPackedArray[etc...}", as a useful way to packing data.

SetDirectory["home/etc.../routes"];
<< airports.mx;
<< routes.mx;

codecoords = 
  Table[ToExpression[airports[[i, 5]]] -> 
    GeoPosition[ToExpression[airports[[i, {-5, -4}]]]], {i, 1, 
    Length[airports]}];
links = Monitor[
   Table[routes[[j, {3, 5}]] /. codecoords, {j, 1, Length[routes]}], 
   ProgressIndicator[j, {1, Length[routes]}]];
linksclean = 
  Developer`ToPackedArray[
   Select[links, 
    Head[#[[1]]] == GeoPosition && Head[#[[2]]] == GeoPosition &]];
With[{locations = RandomChoice[linksclean, 14000]}, 
  GeoGraphics[{{Green, Opacity[0.3], AbsoluteThickness[0.0001], 
     GeoPath[locations, "Geodesic"]}}, GeoRange -> "World", 
   GeoProjection -> Automatic, GeoBackground -> GeoStyling["ReliefMap"], 
   ImageSize -> {1200, 600}]];
Quit[];

Oh dear: 1 in 5 chance that Ebola will spread to the US in September...

and now this: First Ebola case diagnosed in the US

I will post an update of this model very soon.

Cheers, Marco

The substance of this post is great, but I could not permit the following wonderful sentence to go unnoticed.

"Again I had a much more elegant way of doing this, with the advantage that it did not work. "

I'm not sure if you meant advantage or disadvantage, but it is funny and true either way.