# install Pint if necessary

try:
    import pint
except ImportError:
    !pip install pint

# download modsim.py if necessary

from os.path import exists

filename = 'modsim.py'
if not exists(filename):
    from urllib.request import urlretrieve
    url = 'https://raw.githubusercontent.com/AllenDowney/ModSim/main/'
    local, _ = urlretrieve(url+filename, filename)
    print('Downloaded ' + local)

# import functions from modsim

from modsim import *

"""
Explanation: Chapter 8
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
End of explanation
"""


import os

filename = 'World_population_estimates.html'

if not os.path.exists(filename):
    !wget https://raw.githubusercontent.com/AllenDowney/ModSimPy/master/data/World_population_estimates.html

"""
Explanation: In the previous chapter we developed a quadratic model of world
population growth from 1950 to 2016. It is a simple model, but it fits
the data well and the mechanisms it's based on are plausible.
In this chapter we'll use the quadratic model to generate projections of future growth, and compare our results to projections from actual
demographers.
Here's the code that downloads the data.
End of explanation
"""


from pandas import read_html

tables = read_html(filename, header=0, index_col=0, decimal='M')

table2 = tables[2]
table2.columns = ['census', 'prb', 'un', 'maddison', 
                  'hyde', 'tanton', 'biraben', 'mj', 
                  'thomlinson', 'durand', 'clark']

"""
Explanation: And here's the code that reads table2, which contains world populations estimates from the U.S. Census and U.N. DESA, among other organizations.
End of explanation
"""


from modsim import TimeSeries

def run_simulation(system, growth_func):
    results = TimeSeries()
    results[system.t_0] = system.p_0
    
    for t in range(system.t_0, system.t_end):
        growth = growth_func(results[t], t, system)
        results[t+1] = results[t] + growth
        
    return results

def growth_func_quad(pop, t, system):
    return system.alpha * pop + system.beta * pop**2

census = table2.census / 1e9
un = table2.un / 1e9

from modsim import System

t_0 = census.index[0]
p_0 = census[t_0]

system = System(t_0=t_0,
                p_0=p_0,
                t_end=2100)

system.alpha = 25 / 1000
system.beta = -1.8 / 1000

results = run_simulation(system, growth_func_quad)

from modsim import show

show(results.tail())

"""
Explanation: Generating Projections
Now let's run the quadratic model, extending the results until 2100, and see how our projections compare to the professionals'.
Here's the code we'll need from the previous chapter.
End of explanation
"""


from modsim import decorate

results.plot(color='gray', label='model')
decorate(xlabel='Year', 
         ylabel='World population (billion)',
         title='Quadratic Model Projection')

"""
Explanation: And here are the results.
End of explanation
"""


table3 = tables[3]
table3.head()

"""
Explanation: According to the model, population growth will slow gradually after 2020, approaching 12.5 billion by 2100.
I am using the word "projection" deliberately, rather than
"prediction", with the following distinction: "prediction" implies
something like "this is what we should reasonably expect to happen, at
least approximately"; "projection" implies something like "if this
model is actually a good description of what is happening in this
system, and if nothing in the future causes the parameters of the model to change, this is what would happen."
Using "projection" leaves open the possibility that there are important things in the real world that are not captured in the model. It also suggests that, even if the model is good, the parameters we estimate based on the past might be different in the future.
The quadratic model we've been working with is based on the assumption
that population growth is limited by the availability of resources; in
that scenario, as the population approaches carrying capacity, birth
rates fall and death rates rise because resources become scarce.
If that assumption is valid, we might be able to use actual population
growth to estimate carrying capacity, especially if we observe the
transition into the regime where the growth rate starts to fall.
But in the case of world population growth, those conditions don't
apply. Over the last 50 years, the net growth rate has leveled off, but not yet started to fall, so we don't have enough data to make a credible estimate of carrying capacity. And resource limitations are probably not the primary reason growth has slowed. As evidence, consider:


First, the death rate is not increasing; rather, it has declined
    from 1.9% in 1950 to 0.8% now (see http://modsimpy.com/mortality).
    So the decrease in net growth is due entirely to declining birth
    rates.


Second, the relationship between resources and birth rate is the
    opposite of what the model assumes; as nations develop and people
    become more wealthy, birth rates tend to fall.


We should not take too seriously the idea that this model can estimate
carrying capacity. But the predictions of a model can be credible even
if the assumptions of the model are not strictly true. For example,
population growth might behave as if it is resource limited, even if
the actual mechanism is something else.
In fact, demographers who study population growth often use models
similar to ours. In the next section, we'll compare our projections to
theirs.
Projections
From the same page where we got the past population estimates, we'll read table3, which contains predictions for population growth over the next 50-100 years, generated by the U.S. Census, U.N. DESA, and the Population Reference Bureau.
End of explanation
"""


table3.columns = ['census', 'prb', 'un']

"""
Explanation: Some values are NaN, which indicates missing data, because some organizations did not publish projections for some years.
The column names are long strings; for convenience, I'll replace them with abbreviations.
End of explanation
"""


def plot_projections(table):
    """Plot world population projections.
    
    table: DataFrame with columns 'un' and 'census'
    """
    census_proj = table.census.dropna() / 1e9
    un_proj = table.un.dropna() / 1e9
    
    census_proj.plot(style=':', label='US Census')
    un_proj.plot(style='--', label='UN DESA')
    
    decorate(xlabel='Year', 
             ylabel='World population (billion)')

"""
Explanation: The following function plots projections from the U.N. DESA and U.S. Census.  It uses dropna to remove the NaN values from each series before plotting it.
End of explanation
"""


plot_projections(table3)
results.plot(color='gray', label='model')
decorate(title='Quadratic Model Projection')

"""
Explanation: Here are their projections compared to the results of the quadratic model.
End of explanation
"""


diff = census.diff()
diff.head()

"""
Explanation: The U.N. DESA expects the world population to reach 11 billion around 2100, and then level off.
Projections by U.S. Census are a little lower, and they only go until 2050.
Real demographers expect world population to grow more slowly than our model projects, probably because their models are broken down by region and country, where conditions are different, and they take into account expected economic development.
Nevertheless, their projections are qualitatively similar to ours, and
theirs differ from each other almost as much as they differ from ours.
So the results from this model, simple as it is, are not entirely unreasonable.
Summary
You might be interested in this video by Hans Rosling about the demographic changes we expect in this century.
Exercises
Exercise: The net growth rate of world population has been declining for several decades.  That observation suggests one more way to generate more realistic projections, by extrapolating observed changes in growth rate.
To compute past growth rates, we'll use a function called diff, which computes the difference between successive elements in a Series.  For example, here are the changes from one year to the next in census:
End of explanation
"""


alpha = census.diff() / census
alpha.head()

"""
Explanation: The first element is NaN because we don't have the data for 1945, so we can't compute the first difference.
If we divide these differences by the populations, the result is an estimate of the growth rate during each year:
End of explanation
"""


def plot_alpha():
    alpha_census = census.diff() / census
    alpha_census.plot(style='.', label='US Census')

    alpha_un = un.diff() / un
    alpha_un.plot(style='.', label='UN DESA')

    decorate(xlabel='Year', label='Net growth rate')

"""
Explanation: The following function computes and plots the growth rates for the census and un estimates:
End of explanation
"""


plot_alpha()

"""
Explanation: And here's what it looks like.
End of explanation
"""


def alpha_func(t):
    intercept = 0.02
    slope = -0.00021
    return intercept + slope * (t - 1970)

"""
Explanation: Other than a bump around 1990, net growth rate has been declining roughly linearly since 1970.
We can model the decline by fitting a line to this data and extrapolating into the future.
Here's a function that takes a time stamp and computes a line that roughly fits the growth rates since 1970.
End of explanation
"""


from numpy import linspace

t_array = linspace(1960, 2020, 5)
alpha_array = alpha_func(t_array)

"""
Explanation: To see what it looks like, I'll create an array of time stamps from 1960 to 2020 and use alpha_func to compute the corresponding growth rates.
End of explanation
"""


from matplotlib.pyplot import plot

plot_alpha()
plot(t_array, alpha_array, color='gray')

"""
Explanation: Here's what it looks like, compared to the data.
End of explanation
"""


# Solution

t_0 = 1960
t_end = 2100
p_0 = census[t_0]

# Solution

system = System(t_0=t_0, 
                t_end=t_end,
                p_0=p_0,
                alpha_func=alpha_func)

# Solution

def growth_func_alpha(pop, t, system):
    return system.alpha_func(t) * pop

# Solution

growth_func_alpha(p_0, t_0, system)

# Solution

results2 = run_simulation(system, growth_func_alpha);

# Solution

plot_projections(table3)
results2.plot(color='gray', label='model')
decorate(title='Proportional model, linearly decreasing rate')

# Solution

# If the net growth rate continues to decrease linearly,
# world population will peak around 2065 at about 9.8 billion,
# and then start to decline.

# Solution

results.idxmax(), results.max()

"""
Explanation: If you don't like the slope and intercept I chose, feel free to adjust them.
Now, as an exercise, you can use this function to make a projection of world population until 2100.


Create a System object that includes alpha_func as a system variable.


Define a growth function that uses alpha_func to compute the net growth rate at the given time t.


Run a simulation from 1960 to 2100 with your update function, and plot the results.


Compare your projections with those from the US Census and UN.
End of explanation
"""


