Repository: carbonplan/ton-year (language: python)

Issue:
Using an external module to calculate uncertainty
Hey, @jhamman !

There's an old-ish looking repo that has a function to calculate uncertainty for the carbon decay curve that we'd like to use: https://github.com/gschivley/ghgforcing/blob/master/ghgforcing/ghgforcing.py

What do you think is the best way to use this code / create clear attribution? 

A couple of observations: 
- The repo doesn't look active, so it might be nice not to depend on it.  
- Right now, we only want to use the functions related to CO2, not CH4. Would be nice not to bring any mentions of CH4 into the ton-year repo. 

Fork? Copy and paste the functions of interest with attribution? Would love your take!



Discussion / hints:
Great question @freyac - In this case, I think it makes good sense to copy the functions we want into this repo (probably as a clearly separate module), along with the ghgforcing projects License (which is MIT, so all good there). 

Fix patch:
diff --git a/licenses/GHGFORCING b/licenses/GHGFORCING
new file mode 100644
index 0000000..909cee3
--- /dev/null
+++ b/licenses/GHGFORCING
@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2016 Greg Schivley
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/LICENSE b/licenses/LICENSE
similarity index 100%
rename from LICENSE
rename to licenses/LICENSE
diff --git a/requirements.txt b/requirements.txt
index 24ce15a..1c813b6 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,1 +1,3 @@
 numpy
+pandas
+scipy
diff --git a/setup.cfg b/setup.cfg
index 1743d32..d9ec0c6 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -7,7 +7,7 @@ select = B,C,E,F,W,T4,B9
 
 [isort]
 known_first_party=carbonplan
-known_third_party=numpy,pytest,setuptools
+known_third_party=numpy,pandas,pytest,scipy,setuptools
 multi_line_output=3
 include_trailing_comma=True
 force_grid_wrap=0
@@ -16,3 +16,8 @@ line_length=100
 
 [mypy]
 plugins = numpy.typing.mypy_plugin
+
+[mypy-pandas.*]
+ignore_missing_imports = True
+[mypy-scipy.*]
+ignore_missing_imports = True
diff --git a/tonyear/__init__.py b/tonyear/__init__.py
index 5b82742..b05b9c5 100644
--- a/tonyear/__init__.py
+++ b/tonyear/__init__.py
@@ -1,3 +1,4 @@
 # flake8: noqa
 
 from .core import calculate_tonyears, get_baseline_curve, print_benefit_report
+from .ghgforcing import joos_2013, joos_2013_monte_carlo
diff --git a/tonyear/ghgforcing.py b/tonyear/ghgforcing.py
new file mode 100644
index 0000000..fd5e424
--- /dev/null
+++ b/tonyear/ghgforcing.py
@@ -0,0 +1,137 @@
+# Module developed from:
+# https://github.com/gschivley/ghgforcing/blob/master/ghgforcing/ghgforcing.py
+# The MIT License (MIT)
+
+# Copyright (c) 2016 Greg Schivley
+
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+
+import numpy as np
+import pandas as pd
+from scipy.stats import multivariate_normal
+
+
+def joos_2013(t_horizon, **kwargs):
+    """Returns the IRF for CO2 using parameter values from IPCC AR5/Joos et al (2013)
+    Keyword arguments are parameter values.
+
+    Parameters
+    ----------
+    t_horizon : int
+        Length of the time horizon (years)
+
+    Returns
+    -------
+    IRF : np.ndarray
+        IRF curve in the form of an 1D array
+    """
+
+    a0 = kwargs.get('a0', 0.2173)
+    a1 = kwargs.get('a1', 0.224)
+    a2 = kwargs.get('a2', 0.2824)
+    a3 = kwargs.get('a3', 0.2763)
+    tau1 = kwargs.get('tau1', 394.4)
+    tau2 = kwargs.get('tau2', 36.54)
+    tau3 = kwargs.get('tau3', 4.304)
+
+    t = np.arange(t_horizon)
+    IRF = a0 + a1 * np.exp(-t / tau1) + a2 * np.exp(-t / tau2) + a3 * np.exp(-t / tau3)
+    return IRF
+
+
+def joos_2013_monte_carlo(runs: int = 100, t_horizon: int = 1001, **kwargs):
+    """Runs a monte carlo simulation for the Joos_2013 baseline IRF curve
+    using uncertainty parameters for the Joos_2013 curve calculated by
+    Olivie and Peters (2013): https://esd.copernicus.org/articles/4/267/2013/
+
+    Parameters
+    ----------
+    runs : int
+        Number of runs for Monte Carlo simulation. Must be >1.
+
+    t_horizon : int
+        Length of the time horizon over which baseline curve is
+        calculated (years)
+
+    Returns
+    -------
+    summary: pd.DataFrame
+        Dataframe with 'mean', '+sigma', and '-sigma' columns summarizing
+        results of Monte Carlo simulation.
+
+    results: np.array
+        Results from all Monte Carlo runs.
+    """
+
+    if runs <= 1:
+        raise ValueError('number of runs must be >1')
+
+    results = np.zeros((t_horizon, runs))
+
+    # Monte Carlo simulations
+
+    # sigma and x are from Olivie and Peters (2013) Table 5 (J13 values)
+    # They are the covariance and mean arrays for CO2 IRF uncertainty
+    sigma = np.array(
+        [
+            [0.129, -0.058, 0.017, -0.042, -0.004, -0.009],
+            [-0.058, 0.167, -0.109, 0.072, -0.015, 0.003],
+            [0.017, -0.109, 0.148, -0.043, 0.013, -0.013],
+            [-0.042, 0.072, -0.043, 0.090, 0.009, 0.006],
+            [-0.004, -0.015, 0.013, 0.009, 0.082, 0.013],
+            [-0.009, 0.003, -0.013, 0.006, 0.013, 0.046],
+        ]
+    )
+
+    x = np.array([5.479, 2.913, 0.496, 0.181, 0.401, -0.472])
+
+    p_samples = multivariate_normal.rvs(x, sigma, runs)
+    p_df = pd.DataFrame(p_samples, columns=['t1', 't2', 't3', 'b1', 'b2', 'b3'])
+
+    p_exp = np.exp(p_df)
+    a1 = p_exp['b1'] / (1 + p_exp['b1'] + p_exp['b2'] + p_exp['b3'])
+    a2 = p_exp['b2'] / (1 + p_exp['b1'] + p_exp['b2'] + p_exp['b3'])
+    a3 = p_exp['b3'] / (1 + p_exp['b1'] + p_exp['b2'] + p_exp['b3'])
+
+    tau1 = p_exp['t1']
+    tau2 = p_exp['t2']
+    tau3 = p_exp['t3']
+
+    for count in np.arange(runs):
+        co2_kwargs = {
+            'a1': a1[count],
+            'a2': a2[count],
+            'a3': a3[count],
+            'tau1': tau1[count],
+            'tau2': tau2[count],
+            'tau3': tau3[count],
+        }
+
+        irf = joos_2013(t_horizon, **co2_kwargs)
+        results[:, count] = irf
+
+    summary = pd.DataFrame(columns=['mean', '-2sigma', '+2sigma', '5th', '95th'])
+    summary['mean'] = np.mean(results, axis=1)
+    summary['+2sigma'] = summary['mean'] + (1.96 * np.std(results, axis=1))
+    summary['-2sigma'] = summary['mean'] - (1.96 * np.std(results, axis=1))
+    summary['5th'] = np.percentile(results, 5, axis=1)
+    summary['95th'] = np.percentile(results, 95, axis=1)
+
+    return summary, results


Test patch:
diff --git a/tonyear/tests/test_tonyear.py b/tonyear/tests/test_tonyear.py
index d0d140c..648da5c 100644
--- a/tonyear/tests/test_tonyear.py
+++ b/tonyear/tests/test_tonyear.py
@@ -1,7 +1,13 @@
 import numpy as np
 import pytest
 
-from tonyear import calculate_tonyears, get_baseline_curve, print_benefit_report
+from tonyear import (
+    calculate_tonyears,
+    get_baseline_curve,
+    joos_2013,
+    joos_2013_monte_carlo,
+    print_benefit_report,
+)
 
 
 @pytest.mark.parametrize('curve_name', ['joos_2013', 'ipcc_2007', 'ipcc_2000'])
@@ -113,3 +119,27 @@ def test_print_benefit_report():
     }
     ret_val = print_benefit_report(method_dict)
     assert ret_val is None
+
+
+@pytest.mark.parametrize('t_horizon', [1, 100, 1001])
+def test_joos_2013(t_horizon):
+    """
+    Test the the values returned from the Joos 2013 IRF curve creation in the ghgfocing
+    module match the values returned by get_baseline_curve.
+    """
+    curve = joos_2013(t_horizon)
+    assert len(curve) == t_horizon
+    assert np.issubdtype(curve.dtype, np.floating)
+    assert curve.all() == get_baseline_curve('joos_2013', t_horizon=t_horizon).all()
+
+
+def test_joos_2013_monte_carlo_raises_invalid_args():
+    with pytest.raises(ValueError, match='number of runs must be >1'):
+        _ = joos_2013_monte_carlo(1, 1001)
+
+
+@pytest.mark.parametrize('runs', [100, 500])
+@pytest.mark.parametrize('t_horizon', [1, 500, 1001])
+def test_joos_2013_monte_carlo(runs, t_horizon):
+    summary, results = joos_2013_monte_carlo(runs=runs, t_horizon=t_horizon)
+    assert results.shape == (t_horizon, runs)