# How to plot GDPs in this way?

## Actually the current answer is very wonderfull.The bouns just for someone can complete this question.

I found this interesting plot of how the relationship between the GDPs of different countries has changed over time: In Mathematica we can conveniently get historic GDP data like this:

GDP[country_] = CountryData[country, {"GDP", {1973, 2008}}]


How can I recreate the animation above?

• @Dr. belisarius Can you see it now? – yode Mar 24 '16 at 6:34
• Yep. Thanks a lot – Dr. belisarius Mar 24 '16 at 6:35
• – Dr. belisarius Mar 24 '16 at 6:39
• @Dr.belisarius Gratitude to your link.:) – yode Mar 24 '16 at 6:43

This plot is an example of a "force layout". Note that the same thing is happening on two levels. The continents are fighting for space, as are the countries inside the continents. We need to figure out how to implement this fighting for space within a predetermined region.

There are six continents, so first we pick six random points in the unit disk: Then we pretend that each of these points is an electric charge, i.e. they influence each other with a force. If the charge of each point is $q_i$, the charge's position is $r_i$ and $r_i-r_j=r_{ij}$ then $$\mathbf{F} = \frac{q_1 q_2}{|\mathbf{r_{12}}|^2}\mathbf{\hat{r}_{12}}$$ where $\mathbf{\hat{r}}$ is the normalized difference in position. The charges are also constrained to the region such that if they move outside they are immediately replaced in the nearest position in the disk again. Integrating the system we get a new charge configuration: Now imagine again that each charge is a continent. The strength of the charge, $q_i$, is the GDP of the continent. The area of the disk that "belongs" to each continent is its basin of attraction. Finding the basins of attraction we get the following plot: Doing this for each year and using the charge positions of the previous year as the initial positions we can generate the evolution over time. To get the countries inside each continent we can apply the same method to each continent.

## How do we implement something like this, you ask?

Get the data

continents = {
EntityClass["Country", "Africa"],
EntityClass["Country", "Asia"],
EntityClass["Country", "Europe"],
EntityClass["Country", "Oceania"],
EntityClass["Country", "NorthAmerica"],
EntityClass["Country", "SouthAmerica"]
};
SetAttributes[getGDP, Listable]
getGDP[entity_EntityClass, year_] := Total@QuantityMagnitude[
CountryData[#, {"GDP", year}] & /@ EntityList[entity] /. _Missing -> 0
]
getGDP[entity_Entity, year_] := QuantityMagnitude[
CountryData[entity, {"GDP", year}] /. _Missing -> 0
]


Generate initial position for the charges

newPoints[region_, entities_, year_] := Transpose[{
RandomPoint[region, Length@entities],
getGDP[entities, year]/Max[getGDP[entities, year]]
}]


Helper functions for get the position/strength of charges

position[charges_] := First /@ charges;
charge[charges_] := Last /@ charges
updateCharges[charges_, data_] := Transpose[{position[charges], data/Max[data]}]


Compute the net force acting on a charge

This is an implementation of the formula in the introduction.

force[{r1_, q1_}, {r2_, q2_}] := (q1 q2/Norm[r1 - r2]^2) Normalize[r1 - r2]
force[{r1_, q1_}, {r1_, q1_}] := 0
force[charges_] := Total /@ Outer[force, charges, charges, 1]


Implement the constraint

constrain[region_, pts_] := With[{nf = RegionNearest[region]}, nf /@ pts]


Perform Euler integration of charge positions

step[region_, stepSize_][charges_] := Transpose[{
constrain[region, position[charges] + stepSize force[charges]],
charge[charges]
}]
updateChargePositions[region_, stepSize_, steps_, charges_] := Nest[
step[region, stepSize], charges, steps
]


Identify basins of attraction

Basins of attraction are usually computed by dropping sample points into the force field and integrating them to see where they end up (example). Since there is no rotation in the force field and it's very simple, I just drop sample points into the force field and say that each point belongs to the charge which influences it the most. In order to find the boundaries some parts of the force field need to be sampled more carefully. There is an algorithm in PlotRegion which implements adaptive sampling, so I make use of that here.

basinIndex[charges_, particle_] := Last@Ordering[Norm[force[#, {particle, 1}]] & /@ charges]
plotBasins[fullRegion_, charges_] := With[{bounds = RegionBounds[fullRegion]}, Table[RegionPlot[
basinIndex[charges, {x, y}] == i &&
RegionMember[fullRegion, {x, y}],
{x, bounds[[1, 1]], bounds[[1, 2]]},
{y, bounds[[2, 1]], bounds[[2, 2]]}
], {i, Length@continents}]]
getRegions[fullRegion_, charges_] := Cases[Normal[#], p : Polygon[__] :> p, Infinity] & /@ plotBasins[fullRegion, charges]


Make the plot

charges = newPoints[Disk[], continents, 2008];
chargesUpdated = updateChargePositions[Disk[], 10^-2, 3000, charges];
regionsData = getRegions[Disk[], chargesUpdated];
colors = ColorData /@ Range@Length@continents;
edgeform = EdgeForm[{#, Thickness[0.005]}] & /@ colors;
regions = Transpose[{colors, edgeform, regionsData}];
legend = SwatchLegend[ColorData /@ Range@Length@continents, TraditionalForm /@ continents];
Legended[
Graphics[regions, PlotLabel -> "2008"],
legend
]


## Further development

The next step for anyone who is interested is to

1. Write a function that given a year and a set of charges updates the position of those charges to correspond to that year.
2. Apply the function in (1) recursively throughout the years the of interest.
3. Use the charge positions for each year to generate the corresponding plots, and put them together to make an animation.
4. Apply the strategy described here to each continent to divide the continents up into countries.

The machinery more or less exists to perform these for steps. Note that the functions all operate on arbitrary regions, and the continents are specified as regions. This is what makes it possible to recursively apply this method without having to make any major updates. Instead of Disk[] we just need to use RegionUnion @@ continentRegion. RegionUnion is used because each continent is actually made up of many polygons and not just one the way it looks in the plot.

• Interesting.Can you show me How do you get the new charge configuration?? – yode Mar 25 '16 at 3:08
• @yode It's the function updateChargePositions in the section "Perform Euler integration of charge positions" that does that. – C. E. Mar 25 '16 at 3:10
• When data is _Missing ,we can _Missing->(Last year's value) maybe a choice. – yode Mar 25 '16 at 4:01