<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Chapter 7
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License
Step1:  In the previous chapter we developed a population model where net growth during each time step is proportional to the current population. This model seems more realistic than the constant growth model, but it does not fit the data as well.
There are a few things we could try to improve the model
Step3:  And here are the functions from the previous chapter.
Step4:  Quadratic growth
It makes sense that net growth should depend on the current population, but maybe it's not a linear relationship, like this
Step5:  Here's the System object we'll use, initialized with t_0, p_0, and t_end.
Step6:  Now we have to add the parameters alpha and beta .
I chose the following values by trial and error; we'll see better ways to do it later.
Step7:  And here's how we run it
Step8:  And here are the results.
Step9:  The model fits the data well over the whole range, with just a bit of space between them in the 1960s.
Of course, we should expect the quadratic model to fit better than the
constant and proportional models because it has two parameters we can
choose, where the other models have only one. In general, the more
parameters you have to play with, the better you should expect the model
to fit.
But fitting the data is not the only reason to think the quadratic model
might be a good choice. It also makes sense; that is, there is a
legitimate reason to expect the relationship between growth and
population to have this form.
To understand it, let's look at net growth as a function of population.
Net growth
Let's plot the relationship between growth and population in the quadratic model.
I'll use linspace to make an array of 101 populations from 0 to 15 billion.
Step10:  Now I'll use the quadratic model to compute net growth for each population.
Step11:  To plot the growth rate versus population, we can import the plot function from Matplotlib
Step12:  And use it like this.
Step13:  Note that the x-axis is not time, as in the previous figures, but population. We can divide this curve into four regimes of behavior
Step14:  With these parameters, net growth is 0 when the population is about 13.9 billion.
In the context of population modeling, the quadratic model is more
conventionally written like this
Step15:  Now let's see all the ways that can go wrong.
Dysfunction #1
Step16:  This version actually works, but it is not as versatile as it could be.
If there are several System objects, this function can only work with one of them, and only if it is named sys1.
Dysfunction #2
Step17:  In this example, we have a System object named sys1 that gets passed
as an argument to carrying_capacity. But when the function runs, it
ignores the argument and immediately replaces it with a new System
object. As a result, this function always returns the same value, no
matter what argument is passed.
When you write a function, you generally don't know what the values of
the parameters will be. Your job is to write a function that works for
any valid values. If you assign your own values to the parameters, you
defeat the whole purpose of functions.
Dysfunction #3
Step18:  A function that doesn't have a return statement actually returns a special value called None, so in this example the value of pop is None. If you are debugging a program and find that the value of a variable is None when it shouldn't be, a function without a return statement is a likely cause.
Dysfunction #4
Step19:  Exercise
<ASSISTANT_TASK:>
Python Code:
# 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 7
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
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']
un = table2.un / 1e9
census = table2.census / 1e9
Explanation: In the previous chapter we developed a population model where net growth during each time step is proportional to the current population. This model seems more realistic than the constant growth model, but it does not fit the data as well.
There are a few things we could try to improve the model:
Maybe net growth depends on the current population, but the
    relationship is quadratic, not linear.
Maybe the net growth rate varies over time.
In this chapter, we'll explore the first option.
In the exercises, you will have a chance to try the second. 
Here's the code that reads the data.
End of explanation
from modsim import TimeSeries
def run_simulation(system, growth_func):
    Simulate the system using any update function.
    
    system: System object
    growth_func: function that computes the population next year
    
    returns: TimeSeries
    
    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
from modsim import decorate
def plot_estimates():
    census.plot(style=':', label='US Census')
    un.plot(style='--', label='UN DESA')
    decorate(xlabel='Year', 
             ylabel='World population (billion)')
Explanation: And here are the functions from the previous chapter.
End of explanation
def growth_func_quad(pop, t, system):
    return system.alpha * pop + system.beta * pop**2
Explanation: Quadratic growth
It makes sense that net growth should depend on the current population, but maybe it's not a linear relationship, like this:
net_growth = system.alpha * pop
Maybe it's a quadratic relationship, like this:
net_growth = system.alpha * pop + system.beta * pop**2
We can test that conjecture with a new update function:
End of explanation
from modsim import System
t_0 = census.index[0]
p_0 = census[t_0]
t_end = census.index[-1]
system = System(t_0=t_0,
                p_0=p_0,
                t_end=t_end)
Explanation: Here's the System object we'll use, initialized with t_0, p_0, and t_end.
End of explanation
system.alpha = 25 / 1000
system.beta = -1.8 / 1000
Explanation: Now we have to add the parameters alpha and beta .
I chose the following values by trial and error; we'll see better ways to do it later.
End of explanation
results = run_simulation(system, growth_func_quad)
Explanation: And here's how we run it:
End of explanation
results.plot(color='gray', label='model')
plot_estimates()
decorate(title='Quadratic Growth Model')
Explanation: And here are the results.
End of explanation
from numpy import linspace
pop_array = linspace(0, 15, 101)
Explanation: The model fits the data well over the whole range, with just a bit of space between them in the 1960s.
Of course, we should expect the quadratic model to fit better than the
constant and proportional models because it has two parameters we can
choose, where the other models have only one. In general, the more
parameters you have to play with, the better you should expect the model
to fit.
But fitting the data is not the only reason to think the quadratic model
might be a good choice. It also makes sense; that is, there is a
legitimate reason to expect the relationship between growth and
population to have this form.
To understand it, let's look at net growth as a function of population.
Net growth
Let's plot the relationship between growth and population in the quadratic model.
I'll use linspace to make an array of 101 populations from 0 to 15 billion.
End of explanation
growth_array = (system.alpha * pop_array + 
                system.beta * pop_array**2)
Explanation: Now I'll use the quadratic model to compute net growth for each population.
End of explanation
from matplotlib.pyplot import plot
Explanation: To plot the growth rate versus population, we can import the plot function from Matplotlib:
End of explanation
plot(pop_array, growth_array, label='growth')
decorate(xlabel='Population (billions)',
         ylabel='Net growth (billions)',
         title='Growth vs. Population')
Explanation: And use it like this.
End of explanation
-system.alpha / system.beta
Explanation: Note that the x-axis is not time, as in the previous figures, but population. We can divide this curve into four regimes of behavior:
When the population is less than 3-4 billion, net growth is
    proportional to population, as in the proportional model. In this
    regime, the population grows slowly because the population is small.
Between 4 billion and 10 billion, the population grows quickly
    because there are a lot of people.
Above 10 billion, population grows more slowly; this behavior models
    the effect of resource limitations that decrease birth rates or
    increase death rates.
Above 14 billion, resources are so limited that the death rate
    exceeds the birth rate and net growth becomes negative.
Just below 14 billion, there is a point where net growth is 0, which
means that the population does not change. At this point, the birth and death rates are equal, so the population is in equilibrium.
Equilibrium
To find the equilibrium point, we can find the roots, or zeros, of this equation: 
$$\Delta p = \alpha p + \beta p^2$$ 
where $\Delta p$ is net
population growth, $p$ is current population, and $\alpha$ and $\beta$
are the parameters of the model. We can rewrite the right hand side like
this: 
$$\Delta p = p (\alpha + \beta p)$$ 
which is $0$ when $p=0$ or $p=-\alpha/\beta$.
We can use the parameters of the system to compute the second equilibrium point:
End of explanation
def carrying_capacity(system):
    K = -system.alpha / system.beta
    return K
    
sys1 = System(alpha=0.025, beta=-0.0018)
pop = carrying_capacity(sys1)
print(pop)
Explanation: With these parameters, net growth is 0 when the population is about 13.9 billion.
In the context of population modeling, the quadratic model is more
conventionally written like this: 
$$\Delta p = r p (1 - p / K)$$ 
This is the same model; it's just a different way to parameterize it. Given $\alpha$ and $\beta$, we can compute $r=\alpha$ and $K=-\alpha/\beta$.
In this version, it is easier to interpret the parameters: $r$ is the
maximum growth rate, observed when $p$ is small, and $K$ is the
equilibrium point. $K$ is also called the carrying capacity, since
it indicates the maximum population the environment can sustain.
Summary
In this chapter we implemented a quadratic growth model where net growth depends on the current population and the population squared.
This model fits the data well, and we saw one reason why: it is based on the assumption that there is a limit to the number of people the Earth can support.
In the next chapter we'll use the models we have developed to generate
predictions.
But first, I want to warn you about a few things that can go wrong when you write functions.
Dysfunctions
When people learn about functions, there are a few things they often
find confusing. In this section I present and explain some common
problems.
As an example, suppose you want a function that takes a
System object, with variables alpha and beta, and computes the
carrying capacity, -alpha/beta. 
Here's a good solution:
End of explanation
def carrying_capacity():
    K = -sys1.alpha / sys1.beta
    return K
    
sys1 = System(alpha=0.025, beta=-0.0018)
pop = carrying_capacity()
print(pop)
Explanation: Now let's see all the ways that can go wrong.
Dysfunction #1: Not using parameters. In the following version, the
function doesn't take any parameters; when sys1 appears inside the
function, it refers to the object we create outside the function.
End of explanation
# WRONG
def carrying_capacity(system):
    system = System(alpha=0.025, beta=-0.0018)
    K = -system.alpha / system.beta
    return K
    
sys1 = System(alpha=0.03, beta=-0.002)
pop = carrying_capacity(sys1)
print(pop)
Explanation: This version actually works, but it is not as versatile as it could be.
If there are several System objects, this function can only work with one of them, and only if it is named sys1.
Dysfunction #2: Clobbering the parameters. When people first learn
about parameters, they often write functions like this:
End of explanation
# WRONG
def carrying_capacity(system):
    K = -system.alpha / system.beta
    
sys1 = System(alpha=0.025, beta=-0.0018)
pop = carrying_capacity(sys1)
print(pop)
Explanation: In this example, we have a System object named sys1 that gets passed
as an argument to carrying_capacity. But when the function runs, it
ignores the argument and immediately replaces it with a new System
object. As a result, this function always returns the same value, no
matter what argument is passed.
When you write a function, you generally don't know what the values of
the parameters will be. Your job is to write a function that works for
any valid values. If you assign your own values to the parameters, you
defeat the whole purpose of functions.
Dysfunction #3: No return value. Here's a version that computes the
value of K but doesn't return it.
End of explanation
# Solution
system.r = system.alpha
system.K = -system.alpha/system.beta
system.r, system.K
# Solution
def growth_func_quad2(pop, t, system):
    return system.r * pop * (1 - pop / system.K)
# Solution
results2 = run_simulation(system, growth_func_quad2)
results2.plot(color='gray', label='model')
plot_estimates()
decorate(title='Quadratic Growth Model, alternate parameters')
Explanation: A function that doesn't have a return statement actually returns a special value called None, so in this example the value of pop is None. If you are debugging a program and find that the value of a variable is None when it shouldn't be, a function without a return statement is a likely cause.
Dysfunction #4: Ignoring the return value. Finally, here's a version
where the function is correct, but the way it's used is not.
```
def carrying_capacity(system):
    K = -system.alpha / system.beta
    return K
sys1 = System(alpha=0.025, beta=-0.0018)
carrying_capacity(sys1)
print(K)
```
In this example, carrying_capacity runs and returns K, but the
return value doesn't get displayed or assigned to a variable.
If we try to print K, we get a NameError, because K only exists inside the function.
When you call a function that returns a value, you should do something
with the result.
Exercises
Exercise: In a previous section, we saw a different way to parameterize the quadratic model:
$$ \Delta p = r p (1 - p / K) $$
where $r=\alpha$ and $K=-\alpha/\beta$.  
Write a version of growth_func that implements this version of the model.  Test it by computing the values of r and K that correspond to alpha=0.025, beta=-0.0018, and confirm that you get the same results.
End of explanation
# Solution
p0_array = linspace(1, 25, 11)
for p_0 in p0_array:
    system.p_0 = p_0
    results3 = run_simulation(system, growth_func_quad)
    results3.plot(label='_nolegend')
decorate(xlabel='Year',
         ylabel='Population (billions)',
         title='Projections with hypothetical starting populations')
Explanation: Exercise:  What happens if we start with an initial population above the carrying capacity, like 20 billion?  Run the model with initial populations between 1 and 20 billion, and plot the results on the same axes.
Hint: If there are too many labels in the legend, you can plot results like this:
results.plot(label='_nolegend')
End of explanation
<END_TASK>