# 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 10
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

from pandas import read_html

filename = 'World_population_estimates.html'
tables = read_html(filename, header=0, index_col=0, decimal='M')
len(tables)

"""
Explanation: Case studies
This chapter presents case studies where you can practice the tools we have learned so far.
Prehistoric World Population
On the Wikipedia page about world population estimates, the first table contains estimates for prehistoric populations.  The following cells process this table and plot some of the results.
End of explanation
"""


table1 = tables[1]
table1.head()

"""
Explanation: Select tables[1], which is the second table on the page.
End of explanation
"""


table1.columns = ['PRB', 'UN', 'Maddison', 'HYDE', 'Tanton', 
                  'Biraben', 'McEvedy & Jones', 'Thomlinson', 'Durand', 'Clark']

"""
Explanation: Not all researchers provide estimates for the same dates.
Again, we'll replace the long column names with more convenient abbreviations.
End of explanation
"""


for col in table1.columns:
    table1[col] = pd.to_numeric(table1[col], errors='coerce')

"""
Explanation: Some of the estimates are in a form Pandas doesn't recognize as numbers, but we can coerce them to be numeric.
End of explanation
"""


table1.plot()
decorate(xlim=[-10000, 2000], xlabel='Year', 
         ylabel='World population (millions)',
         title='Prehistoric population estimates')
plt.legend(fontsize='small');

"""
Explanation: Here are the results.  Notice that we are working in millions now, not billions.
End of explanation
"""


table1.plot()
decorate(xlim=[0, 2000], xlabel='Year', 
         ylabel='World population (millions)',
         title='CE population estimates')

"""
Explanation: We can use xlim to zoom in on everything after Year 0.
End of explanation
"""


# Solution

# The function I found that best matches the data has the form
# a + b / (c - x)

# This function is hard to explain physically; that is, it doesn't
# correspond to a growth model that makes sense in terms of human 
# behavior.

# And it implies that the population goes to infinity in 2040.

from numpy import linspace
from matplotlib.pyplot import plot

xs = linspace(100, 1950)
ys = 110 + 200000 / (2040 - xs)
table1.plot()
plot(xs, ys, color='gray', alpha=0.5, label='model')

decorate(xlim=[0, 2000], xlabel='Year', 
         ylabel='World population (millions)',
         title='CE population estimates')

# Solution

# And it doesn't do a particularly good job of predicting
# actual growth from 1940 to the present.

# TODO: load these series from table2
#plot(census, ':', label='US Census')
#plot(un, '--', label='UN DESA')

xs = linspace(1940, 2020)
ys = 110 + 200000 / (2040 - xs)
plot(xs, ys/1000, color='gray', label='model')

decorate(xlim=[1950, 2016], xlabel='Year', 
         ylabel='World population (billions)',
         title='Prehistoric population estimates')

"""
Explanation: See if you can find a model that fits these data well from Year 0 to 1950.
How well does your best model predict actual population growth from 1950 to the present?
End of explanation
"""


