<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Trees
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License
Step1:  Modeling tree growth
This case study is based on "Height-Age Curves for Planted Stands of Douglas Fir, with Adjustments for Density", a working paper by Flewelling, Collier, Gonyea, Marshall, and Turnblom.
It provides "site index curves", which are curves that show the expected height of the tallest tree in a stand of Douglas firs as a function of age, for a stand where the trees are the same age.
Depending on the quality of the site, the trees might grow more quickly or slowly.  So each curve is identified by a "site index" that indicates the quality of the site.
I'll start with some of the data from their Table 1.  Here's the sequence of ages.
Step2:  And here's the series of heights for a site with index 45, indicating that height at 30 years is 45 feet.
Step3:  Here's the series for site index 65.
Step4:  And for site index 85.
Step5:  Here's what the curves look like
Step6:  For my examples I'll work with the SI 65 data; as an exercise, you can run the notebook again with either of the other curves.
Step7:  Model 1
As a starting place, let's assume that the ability of the tree to gain mass is limited by the area it exposes to sunlight, and that the growth rate (in mass) is proportional to that area.  In that case we can write
Step9:  And here's an update function that takes the current height as a parameter and returns the height during the next time step.
Step10:  I'll test the update function with the initial conditions.
Step12:  Here's our usual version of run_simulation.
Step13:  And here's how we run it.
Step14:  Here's what the results look like
Step15:  The model converges to a straight line.
I chose the value of alpha to fit the data as well as I could, but it is clear that the data have curvature that's not captured by the model.
Here are the errors, that is, the differences between the model and the data.
Step16:  And here's the mean absolute error.
Step17:  This model might explain why the height of a tree grows roughly linearly
Step19:  I'll wrap the code from the previous section is a function that takes the parameters as inputs and makes a System object.
Step20:  Here's how we use it.
Step21:  With different values for the parameters, we get curves with different behavior.  Here are a few that I chose by hand.
Step23:  To find the parameters that best fit the data, I'll use leastsq.
We need an error function that takes parameters and returns errors
Step24:  Here's how we use it
Step25:  Now we can pass error_func to leastsq, which finds the parameters that minimize the squares of the errors.
Step26:  Using the best parameters we found, we can run the model and plot the results.
Step27:  The mean absolute error is better than for Model 1, but that doesn't mean much.  The model still doesn't fit the data well.
Step28:  And the estimated fractal dimension is 3.11, which doesn't seem likely.
Let's try one more thing.
Model 3
Models 1 and 2 imply that trees can grow forever, but we know that's not true.  As trees get taller, it gets harder for them to move water and nutrients against the force of gravity, and their growth slows.
We can model this effect by adding a term to the model similar to what we saw in the logistic model of population growth.  Instead of assuming
Step30:  Here's an updated version of make_system
Step31:  Here's the new System object.
Step33:  And here's the new update function.
Step34:  As always, we'll test the update function with the initial conditions.
Step35:  And we'll test the error function with the new update function.
Step36:  Now let's search for the best parameters.
Step37:  With these parameters, we can fit the data much better.
Step38:  And the mean absolute error is substantially smaller.
Step39:  The estimated fractal dimension is about 2.6, which is plausible; it suggests that if you double the height of the tree, the mass grows by a factor of $2^{2.6}$
<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: Trees
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
End of explanation
years = [2, 3, 4, 5, 6, 8, 10, 15, 20, 25, 30,
         35, 40, 45, 50, 55, 60, 65, 70]
Explanation: Modeling tree growth
This case study is based on "Height-Age Curves for Planted Stands of Douglas Fir, with Adjustments for Density", a working paper by Flewelling, Collier, Gonyea, Marshall, and Turnblom.
It provides "site index curves", which are curves that show the expected height of the tallest tree in a stand of Douglas firs as a function of age, for a stand where the trees are the same age.
Depending on the quality of the site, the trees might grow more quickly or slowly.  So each curve is identified by a "site index" that indicates the quality of the site.
I'll start with some of the data from their Table 1.  Here's the sequence of ages.
End of explanation
site45 = TimeSeries([1.4, 1.49, 1.75, 2.18, 2.78, 4.45, 6.74,
                    14.86, 25.39, 35.60, 45.00, 53.65, 61.60,
                    68.92, 75.66, 81.85, 87.56, 92.8, 97.63],
                    index=years)
Explanation: And here's the series of heights for a site with index 45, indicating that height at 30 years is 45 feet.
End of explanation
site65 = TimeSeries([1.4, 1.56, 2.01, 2.76, 3.79, 6.64, 10.44, 
                    23.26, 37.65, 51.66, 65.00, 77.50, 89.07, 
                    99.66, 109.28, 117.96, 125.74, 132.68, 138.84],
                    index=years)
Explanation: Here's the series for site index 65.
End of explanation
site85 = TimeSeries([1.4, 1.8, 2.71, 4.09, 5.92, 10.73, 16.81, 
                 34.03, 51.26, 68.54, 85, 100.34, 114.33,
                 126.91, 138.06, 147.86, 156.39, 163.76, 170.10],
               index=years)
Explanation: And for site index 85.
End of explanation
site85.plot(label='SI 85')
site65.plot(label='SI 65')
site45.plot(label='SI 45')
decorate(xlabel='Time (years)',
         ylabel='Height (feet)')
Explanation: Here's what the curves look like:
End of explanation
data = site65
Explanation: For my examples I'll work with the SI 65 data; as an exercise, you can run the notebook again with either of the other curves.
End of explanation
alpha = 7
dim = 3
t_0 = data.index[0]
h_0 = data[t_0]
t_end = data.index[-1]
system = System(alpha=alpha, 
                dim=dim, 
                h_0=h_0, 
                t_0=t_0, 
                t_end=t_end)
Explanation: Model 1
As a starting place, let's assume that the ability of the tree to gain mass is limited by the area it exposes to sunlight, and that the growth rate (in mass) is proportional to that area.  In that case we can write:
$$ m_{n+1} = m_n + \alpha A$$
where $m_n$ is the mass of the at time step $n$, $A$ is the area exposed to sunlight, and $\alpha$ is an unknown growth parameter.
To get from $m$ to $A$, I'll make the additional assumption that mass is proportional to height raised to an unknown power:
$$ m = \beta h^D $$
where $h$ is height, $\beta$ is an unknown constant of proportionality, and $D$ is the dimension that relates height and mass. We'll start by assuming $D=3$, but we'll revisit that assumption.
Finally, we'll assume that area is proportional to height squared:
$$ A = \gamma h^2$$
I'll specify height in feet, and choose units for mass and area so that $\beta=1$ and $\gamma=1$.
Putting all that together, we can write a difference equation for height:
$$ h_{n+1}^D = h_n^D + \alpha h_n^2 $$
Now let's simulate this system.  Here's a system object with the parameters and initial conditions.
End of explanation
def update(height, t, system):
    Update height based on geometric model.
    
    height: current height in feet
    t: what year it is
    system: system object with model parameters
    
    area = height**2
    mass = height**system.dim
    mass += system.alpha * area
    return mass**(1/system.dim)
Explanation: And here's an update function that takes the current height as a parameter and returns the height during the next time step.
End of explanation
update(h_0, t_0, system)
Explanation: I'll test the update function with the initial conditions.
End of explanation
def run_simulation(system, update_func):
    Simulate the system using any update function.
    
    system: System object
    update_func: function that computes the population next year
    
    returns: TimeSeries
    
    results = TimeSeries()
    results[system.t_0] = system.h_0
    
    for t in linrange(system.t_0, system.t_end-1):
        results[t+1] = update_func(results[t], t, system)
        
    return results
Explanation: Here's our usual version of run_simulation.
End of explanation
results = run_simulation(system, update)
results.tail()
Explanation: And here's how we run it.
End of explanation
def plot_results(results, data):
    results.plot(style=':', label='model', color='gray')
    data.plot(label='data')
    decorate(xlabel='Time (years)',
             ylabel='Height (feet)')
    
plot_results(results, data)
Explanation: Here's what the results look like:
End of explanation
errors = results - data
errors.dropna()
Explanation: The model converges to a straight line.
I chose the value of alpha to fit the data as well as I could, but it is clear that the data have curvature that's not captured by the model.
Here are the errors, that is, the differences between the model and the data.
End of explanation
def mean_abs_error(results, data):
    return (results-data).abs().mean()
mean_abs_error(results, data)
Explanation: And here's the mean absolute error.
End of explanation
alpha = 7
dim = 2.5
Explanation: This model might explain why the height of a tree grows roughly linearly:
If area is proportional to $h^2$ and mass is proportional to $h^3$, and
Change in mass is proportional to area, and
Height grows linearly, then
Area grows in proportion to $h^2$, and
Mass grows in proportion to $h^3$.
If the goal is to explain (approximate) linear growth, we might stop there.  But this model does not fit the data particularly well, and it implies that trees could keep growing forever.
So we might want to do better.
Model 2
As a second attempt, let's suppose that we don't know $D$.  In fact, we don't, because trees are not like simple solids; they are more like fractals, which have fractal dimension.
I would expect the fractal dimension of a tree to be between 2 and 3, so I'll guess 2.5.
End of explanation
def make_system(params, data):
    Makes a System object.
    
    params: sequence of alpha, dim
    data: Series
    
    returns: System object
    
    alpha, dim = params
    
    t_0 = data.index[0]
    t_end = data.index[-1]
    h_0 = data[t_0]
    return System(alpha=alpha, dim=dim, 
                  h_0=h_0, t_0=t_0, t_end=t_end)
Explanation: I'll wrap the code from the previous section is a function that takes the parameters as inputs and makes a System object.
End of explanation
params = alpha, dim
system = make_system(params, data)
Explanation: Here's how we use it.
End of explanation
def run_and_plot(alpha, dim, data):
    params = alpha, dim
    system = make_system(params, data)
    results = run_simulation(system, update)
    results.plot(style=':', color='gray', label='_nolegend')
run_and_plot(0.145, 2, data)
run_and_plot(0.58, 2.4, data)
run_and_plot(2.8, 2.8, data)
run_and_plot(6.6, 3, data)
run_and_plot(15.5, 3.2, data)
run_and_plot(38, 3.4, data)
data.plot(label='data')
decorate(xlabel='Time (years)',
             ylabel='Height (feet)')
Explanation: With different values for the parameters, we get curves with different behavior.  Here are a few that I chose by hand.
End of explanation
def error_func(params, data, update_func):
    Runs the model and returns errors.
    
    params: sequence of alpha, dim
    data: Series
    update_func: function object
    
    returns: Series of errors
    
    print(params)
    system = make_system(params, data)
    results = run_simulation(system, update_func)
    return (results - data).dropna()
Explanation: To find the parameters that best fit the data, I'll use leastsq.
We need an error function that takes parameters and returns errors:
End of explanation
errors = error_func(params, data, update)
Explanation: Here's how we use it:
End of explanation
best_params, details = leastsq(error_func, params, data, update)
print(details.success)
Explanation: Now we can pass error_func to leastsq, which finds the parameters that minimize the squares of the errors.
End of explanation
system = make_system(best_params, data)
results = run_simulation(system, update)
plot_results(results, data)
Explanation: Using the best parameters we found, we can run the model and plot the results.
End of explanation
mean_abs_error(results, data)
Explanation: The mean absolute error is better than for Model 1, but that doesn't mean much.  The model still doesn't fit the data well.
End of explanation
alpha = 2.0
dim = 2.5
K = 150
params = [alpha, dim, K]
Explanation: And the estimated fractal dimension is 3.11, which doesn't seem likely.
Let's try one more thing.
Model 3
Models 1 and 2 imply that trees can grow forever, but we know that's not true.  As trees get taller, it gets harder for them to move water and nutrients against the force of gravity, and their growth slows.
We can model this effect by adding a term to the model similar to what we saw in the logistic model of population growth.  Instead of assuming:
$ m_{n+1} = m_n + \alpha A $ 
Let's assume
$ m_{n+1} = m_n + \alpha A (1 - h / K) $
where $K$ is similar to the carrying capacity of the logistic model.  As $h$ approaches $K$, the factor $(1 - h/K)$ goes to 0, causing growth to level off.
Here's what the implementation of this model looks like:
End of explanation
def make_system(params, data):
    Makes a System object.
    
    params: sequence of alpha, dim, K
    data: Series
    
    returns: System object
    
    alpha, dim, K = params
    
    t_0 = data.index[0]
    t_end = data.index[-1]
    h_0 = data[t_0]
    return System(alpha=alpha, dim=dim, K=K, 
                  h_0=h_0, t_0=t_0, t_end=t_end)
Explanation: Here's an updated version of make_system
End of explanation
system = make_system(params, data)
Explanation: Here's the new System object.
End of explanation
def update3(height, t, system):
    Update height based on geometric model with growth limiting term.
    
    height: current height in feet
    t: what year it is
    system: system object with model parameters
    
    area = height**2
    mass = height**system.dim
    mass += system.alpha * area * (1 - height/system.K)
    return mass**(1/system.dim)
Explanation: And here's the new update function.
End of explanation
update3(h_0, t_0, system)
Explanation: As always, we'll test the update function with the initial conditions.
End of explanation
error_func(params, data, update3)
Explanation: And we'll test the error function with the new update function.
End of explanation
best_params, details = leastsq(error_func, params, data, update3)
details.success
Explanation: Now let's search for the best parameters.
End of explanation
system = make_system(best_params, data)
results = run_simulation(system, update3)
plot_results(results, data)
Explanation: With these parameters, we can fit the data much better.
End of explanation
mean_abs_error(results, data)
Explanation: And the mean absolute error is substantially smaller.
End of explanation
2**2.6
Explanation: The estimated fractal dimension is about 2.6, which is plausible; it suggests that if you double the height of the tree, the mass grows by a factor of $2^{2.6}$
End of explanation
<END_TASK>