4
$\begingroup$

I was looking at the Mathematica code given at http://community.wolfram.com/groups/-/m/t/326240

where example code is given to simulate the spreading of a supposed Ebola outbreak in Nigeria:

CountryData["Nigeria", "Population"]
Graphics[CountryData["Nigeria", "Polygon"]]

names = CityData[{All, "Nigeria"}];
citypop = 
  Table[CityData[names[[i]], "Population"], {i, 1, Length[names]}];
citycoords = 
  Table[CityData[names[[i]], "Coordinates"], {i, 1, Length[names]}];

Needs["ComputationalGeometry`"]
dtri = DelaunayTriangulation[citycoords]; list = {}; Table[
 Do[AppendTo[list, {i, dtri[[All, 2]][[i, j]]}], {j, 1, 
   Length[dtri[[All, 2]][[i, All]]]}], {i, 1, Length[dtri]}];
coupling = Table[0, {i, 1, Length[names]}, {j, 1, Length[names]}];
For[i = 1, i < Length[list] + 1, i++, 
 coupling[[list[[i]][[1]], list[[i]][[2]]]] = 1;]

coulinginterm = coupling;
coupling = SparseArray[coulinginterm];

Graphics[Join[
  Table[Circle[citycoords[[i]], 0.02], {i, 1, Length[names]}], 
  DeleteCases[
   Flatten[Table[
     If[coupling[[i, j]] == 1, 
      Line[{citycoords[[i]], citycoords[[j]]}]], {i, 1, 
      Length[names]}, {j, 1, i}], 1], Null]]]

\[Rho] = 0.2; \[Lambda] = 0.1; Mcities = Length[names]; \[Mu] = 0.05;

For[i = 1, i <= Mcities, i++, sus[1, i] = 1.; inf[1, i] = 0.0;
 rec[1, i] = 0.;]

sus[1, 1] = 0.95; inf[1, 1] = 0.05;
rec[1, 1] = 0.0;

meanNN = coupling.citypop;

sumind = Table[
   Take[Flatten[ArrayRules[coupling[[k, All]]][[All, 1]]], 
    Length[ArrayRules[coupling[[k, All]]]] - 1], {k, 1, Mcities}];

sus[i_, j_] := 
  sus[i, j] = (1 - \[Mu]) (sus[i - 1, 
        j] - \[Rho] sus[i - 1, j] inf[i - 1, j]) + \[Mu] Total[
       Table[sus[i - 1, sumind[[j, u]]]*citypop[[sumind[[j, u]]]], {u,
          1, Length[sumind[[j]]]}]]/meanNN[[j]];
inf[i_, j_] := 
  inf[i, j] = (1 - \[Mu]) (inf[i - 1, 
        j] + \[Rho] sus[i - 1, j] inf[i - 1, j] - \[Lambda] inf[i - 1,
          j]) + \[Mu] Total[
       Table[inf[i - 1, sumind[[j, u]]]*citypop[[sumind[[j, u]]]], {u,
          1, Length[sumind[[j]]]}]]/meanNN[[j]];
rec[i_, j_] := 
  rec[i, j] = (1 - \[Mu]) (rec[i - 1, 
        j] + \[Lambda] inf[i - 1, j]) + \[Mu] Total[
       Table[rec[i - 1, sumind[[j, u]]]*citypop[[sumind[[j, u]]]], {u,
          1, Length[sumind[[j]]]}]]/meanNN[[j]];

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

poly = Graphics[
   Polygon[Flatten[CountryData["Nigeria", "Coordinates"], 1]], 
   ImagePadding -> None];
Animate[ImageSubtract[
  Graphics[ListDensityPlot[
    Join[{{4, 3.25, 0}, {4, 14, 0}, {14, 3.25, 0}, {14, 14, 0}}, 
     Table[{citycoords[[k]][[1]], citycoords[[k]][[2]], 
       1. - tcourse[[t, k, 1]]}, {k, 1, Mcities}]], 
    InterpolationOrder -> 3, ColorFunction -> "Rainbow", 
    PlotRange -> All, Frame -> False, PlotRangePadding -> None]], 
  poly], {t, 1, 500, 1}, DefaultDuration -> 20.]

enter image description here

In Mathematica 9 this runs just fine, and the line

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

just takes 29 seconds to complete the calculation of the system of recurrence equations.

If I try the same in Mathematica 10 it takes ages to complete though - in fact I didn't have the patience to wait for it to complete... Does anyone know what could be wrong in Mathematica 10, or what needs changing to make it run at the same speed as in Mathematica 9?

$\endgroup$
  • 2
    $\begingroup$ I think the problem is most probably not with solving the recurrence equations as such but the fact that CityData and CountryData do return Entity and Quantity objects since version 10. While it is possible to use these in many calculations they of course are slowing down things remarkably. So I would suggest to get rid of those using EntityValue and QuantityValue. I think there is also a global option to change the behavior of the data functions to what it was in version 9, to be found on this site... $\endgroup$ – Albert Retey Dec 9 '14 at 11:19
  • 1
    $\begingroup$ Many thanks - SetSystemOptions["DataOptions" -> "ReturnQuantities" -> False] indeed seems to solve the issue using a global option. How would I have to use EntityValue to get names = CityData[{All, country}]; citypop = Table[CityData[names[[i]], "Population"], {i, 1, Length[names]}]; citycoords = Table[CityData[names[[i]], "Coordinates"], {i, 1, Length[names]}]; working though without setting global options? $\endgroup$ – Tom Wenseleers Dec 9 '14 at 11:43
  • $\begingroup$ Please see my answer for some examples about how to get rid of Quantity and Entity in your data without changing global options... $\endgroup$ – Albert Retey Dec 9 '14 at 12:11
  • 1
    $\begingroup$ Possible duplicate $\endgroup$ – bobthechemist Dec 9 '14 at 13:05
6
$\begingroup$

It looks like the real problem is not with solving the recurrence relation as such but that the CityData and CountryData functions do return some of their results differently in version 10. Many of these return values are now not simple numbers or strings anymore but rather Entity or Quantity "objects". While it is in many cases possible to use these as input to other functions and e.g. use Quantity expressions in many calculations such calculations will become very slow when using Quantity objects instead of just number. The cure then would be to get rid of the Quantity end Entity objects to get speed that is comparable to version 9.

To my understanding the recommended way to get rid of Entity would be to use something like EntityValue[#, "Name"] & /@ CityData[{All, "Nigeria"}] but it turns out that this is really slow and even caused some network timeouts on my computer. Alternatively you can just extract the city "definitions" as lists of strings like this: CityData[{All, "Nigeria"}][[All, 2]] or -- if necessary -- the names (which might not be unique in general) with CityData[{All, "Nigeria"}][[All, 2, 1]]. You can then extract the population numbers which you have to convert to numbers with QuantityMagnitude and the coordinates as follows:

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

While conceptionally the new strategy to use quantities and entities for these kind of data looks very reasonable the current implementation seems to need improvement to be useful in practice, but this might be not too surprising considering this is all brand new stuff. Not only are some very simple operations very slow, but there are also some inconsistencies, as e.g. it is not clear why CityData["Lagos","Coordinates"] returns just a list of numbers while CountryData["Nigeria", "Polygon"] returns a polygon of GeoPositions. For the moment the best way to handle this seems to be to convert such data to more efficient representations and do that with making use of the internal representations instead of official documented functions where necessary. I'm afraid that you might have to adopt these parts of the code for the next versions again...

EDIT: the best way to extract the city definitions as a list of strings seems to be the use of CanonicalName:

cities = CanonicalName[CityData[{All, "Nigeria"}]];

it is as fast as the part variant, is listable and it should not depend on implementation details and thus probably also work well in future versions. It also should work with any other kind of Entity expressions, where it depends on the entity class what exactly the return value is (string vs. list of string etc.)...

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Thx millions for this - that works great!! I also added names = cities[[All, 1]] $\endgroup$ – Tom Wenseleers Dec 9 '14 at 13:00

Not the answer you're looking for? Browse other questions tagged or ask your own question.