Repository: xarray-contrib/xskillscore

Issue:
Add MAPE in See Also and in testing from sklearn
`sklearn` has MAPE now in their master branch (https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/metrics/_regression.py#L197). Once released, we can test against this and put it in "See Also" for MAPE.

Discussion:
raybellwaves: believe the python 3.6 issues in the CI are outside this PR. the doc yaml is set at 3.6 which is probably why it's also failing.
raybellwaves: > can you also add this epsilon to other deterministic metrics? or should this be another PR?

I vote leave for another PR. Could create an issue and list which metrics can be updated.

> tests are failing

Yes deffo something weird with numba and py3.6
```
/usr/share/miniconda/envs/xskillscore-minimum-tests/site-packages/numba/core/config.py:20: in <module>
    MACHINE_BITS = tuple.__itemsize__ * 8
E   AttributeError: type object 'tuple' has no attribute '__itemsize__'
=============================== warnings summary ==========================
```

I may move the doc yaml file to 3.* or 3.8 as it's currently as 3.6

> maybe you could also add a test for the epsilon case

Maybe not be a bad idea. I Could test with a DataArray with some 0's sprinkled in. What would be the best way to do this? create an `o_with_zeros` in the conftest?
aaronspring: > Maybe not be a bad idea. I Could test with a DataArray with some 0's sprinkled in. What would be the best way to do this? create an o_with_zeros in the conftest?

yes
raybellwaves: @aaronspring  What's the syntax for parametrizing a variable from the conftest?

e.g. trying

```
@pytest.mark.parametrize("a_1d", [a_1d, a_1d_with_zeros])
def test_xs_same_as_skl_mape(a_1d, b_1d):
    actual = mape(a_1d, b_1d, "time")
    expected = mean_absolute_percentage_error(a_1d, b_1d)
    assert np.allclose(actual, expected)
```

but says a_1d is not defined. It picks up a_1d in
```
def test_xs_same_as_skl_mape(a_1d, b_1d):
    actual = mape(a_1d, b_1d, "time")
    expected = mean_absolute_percentage_error(a_1d, b_1d)
    assert np.allclose(actual, expected)
```

I guess in the parametrize line it doesn't know `a_1d` and `a_1d_with_zeros` come from the conftest
raybellwaves: I also changed some other things to possible help with the CI. Where we had a GHA testing again 3.6 I moved it to 3.8. I also updated the doc yaml to 3.* from 3.6.
raybellwaves: I rewrote the test so I think i'm good there
aaronspring: Probably fix this numba thing in a separate PR as it persists with 3.6.

maybe just taking a more recent version works, although I didn’t find the cause in https://github.com/numba/numba/blob/master/CHANGE_LOG
aaronspring: `MACHINE_BITS = tuple.__itemsize__ * 8`
https://github.com/numba/numba/blob/ba0a6c8e7a7e33ad812bc644b29d2dfd8f3f4f0d/numba/core/config.py#L20

cannot reproduce locally. `tuple` is from python, so there seems to be some problem in GHA? 
raybellwaves: > cannot reproduce locally. tuple is from python, so there seems to be some problem in GHA?

Yeah it's probably the architecture GHA runs on e.g. linux. I believe there is a place where you can report it for the GHA team to look at it.

We can either merge this knowing that the py3.6 CI is bogus or remove py3.6 from the test suite.
raybellwaves: @aaronspring if you don't mind, for the future can you do your test commits on a branch off of mine rather than this PR
aaronspring: what do we do with that failing py36 testing?
raybellwaves: > what do we do with that failing py36 testing?

If you are happy to we can merge this.

I'll create a separate issue for the py3.6 CI failing.
aaronspring: Just remove my small try in the squash merge
raybellwaves: > Just remove my small try in the squash merge

I think it's ok. I don't really care about the commit history. As long as the files changed are desired.

Patch:
diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml
index 1be599e..8c0f2fe 100644
--- a/.github/workflows/lint.yml
+++ b/.github/workflows/lint.yml
@@ -13,6 +13,6 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v2
         with:
-          python-version: 3.6
+          python-version: 3.8
       - name: Lint via pre-commit checks
         uses: pre-commit/action@v2.0.0
diff --git a/.github/workflows/xskillscore_installs.yml b/.github/workflows/xskillscore_installs.yml
index 68dfc61..4dc3d9d 100644
--- a/.github/workflows/xskillscore_installs.yml
+++ b/.github/workflows/xskillscore_installs.yml
@@ -17,7 +17,7 @@ jobs:
     - name: Setup python
       uses: actions/setup-python@v2
       with:
-        python-version: 3.6
+        python-version: 3.8
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index e0b8757..211bb5b 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -34,6 +34,8 @@ Internal Changes
 - Added Python 3.7 and Python 3.8 to the CI. Use the latest version of Python 3
   for development. (:issue:`21`, :pr:`189`). `Aaron Spring`_
 - Lint with the latest black. (:issue:`179`, :pr:`191`). `Ray Bell`_
+- Update mape algorithm from scikit-learn v0.24.0 and test against it.
+  (:issue:`160`, :pr:`230`) `Ray Bell`_
 
 
 xskillscore v0.0.18 (2020-09-23)
diff --git a/ci/doc.yml b/ci/doc.yml
index afbdf91..ee45765 100644
--- a/ci/doc.yml
+++ b/ci/doc.yml
@@ -2,7 +2,7 @@ name: xskillscore-docs
 channels:
   - conda-forge
 dependencies:
-  - python=3.6
+  - python=3.8
   - bottleneck
   - doc8
   - importlib_metadata
diff --git a/xskillscore/core/deterministic.py b/xskillscore/core/deterministic.py
index aafe46a..a8958eb 100644
--- a/xskillscore/core/deterministic.py
+++ b/xskillscore/core/deterministic.py
@@ -1046,12 +1046,13 @@ def mape(a, b, dim=None, weights=None, skipna=False, keep_attrs=False):
     .. math::
         \\mathrm{MAPE} = \\frac{1}{n} \\sum_{i=1}^{n}
                          \\frac{\\vert a_{i} - b_{i} \\vert}
-                               {\\vert a_{i} \\vert}
+                               {max(\epsilon, \\vert a_{i} \\vert)}
 
     .. note::
         The percent error is calculated in reference to ``a``. Percent
         error is reported as decimal percent. I.e., a value of 1 is
-        100%.
+        100%. :math:`\epsilon` is an arbitrary small yet strictly positive
+        number to avoid undefined results when ``a`` is zero.
 
     Parameters
     ----------
@@ -1078,6 +1079,10 @@ def mape(a, b, dim=None, weights=None, skipna=False, keep_attrs=False):
     xarray.Dataset or xarray.DataArray
         Mean Absolute Percentage Error.
 
+    See Also
+    --------
+    sklearn.metrics.mean_absolute_percentage_error
+
     References
     ----------
     https://en.wikipedia.org/wiki/Mean_absolute_percentage_error
diff --git a/xskillscore/core/np_deterministic.py b/xskillscore/core/np_deterministic.py
index e6fb19a..ae736ac 100644
--- a/xskillscore/core/np_deterministic.py
+++ b/xskillscore/core/np_deterministic.py
@@ -652,7 +652,7 @@ def _mape(a, b, weights, axis, skipna):
     .. math::
         \\mathrm{MAPE} = \\frac{1}{n} \\sum_{i=1}^{n}
                          \\frac{\\vert a_{i} - b_{i} \\vert}
-                               {\\vert a_{i} \\vert}
+                               {max(\epsilon, \\vert a_{i} \\vert)}
 
     Parameters
     ----------
@@ -679,6 +679,13 @@ def _mape(a, b, weights, axis, skipna):
     Percent error is reported as decimal percent. I.e., a value of
     1 is 100%.
 
+    \epsilon is an arbitrary small yet strictly positive number to avoid
+    undefined results when ``a`` is zero.
+
+    See Also
+    --------
+    sklearn.metrics.mean_absolute_percentage_error
+
     References
     ----------
     https://en.wikipedia.org/wiki/Mean_absolute_percentage_error
@@ -687,8 +694,8 @@ def _mape(a, b, weights, axis, skipna):
     if skipna:
         a, b, weights = _match_nans(a, b, weights)
     weights = _check_weights(weights)
-    # replace divided by 0 with nan
-    mape = np.absolute(a - b) / np.absolute(np.where(a != 0, a, np.nan))
+    epsilon = np.finfo(np.float64).eps
+    mape = np.absolute(a - b) / np.maximum(np.absolute(a), epsilon)
     if weights is not None:
         return sumfunc(mape * weights, axis=axis) / sumfunc(weights, axis=axis)
     else:


Test patch:
diff --git a/.github/workflows/xskillscore_testing.yml b/.github/workflows/xskillscore_testing.yml
index 4ed6d27..0508ea9 100644
--- a/.github/workflows/xskillscore_testing.yml
+++ b/.github/workflows/xskillscore_testing.yml
@@ -60,7 +60,7 @@ jobs:
           channels: conda-forge
           mamba-version: '*'
           activate-environment: xskillscore-docs-notebooks
-          python-version: 3.6
+          python-version: 3.8
       - name: Set up conda environment
         run: |
           mamba env update -f ci/docs_notebooks.yml
diff --git a/xskillscore/tests/conftest.py b/xskillscore/tests/conftest.py
index 059f3e2..9998ae9 100644
--- a/xskillscore/tests/conftest.py
+++ b/xskillscore/tests/conftest.py
@@ -68,6 +68,13 @@ def b_nan(b):
     return b.copy().where(b < 0.5)
 
 
+# with zeros
+@pytest.fixture
+def a_with_zeros(a):
+    """Zeros"""
+    return a.copy().where(a < 0.5, 0)
+
+
 # dask
 @pytest.fixture
 def a_dask(a):
@@ -116,6 +123,12 @@ def b_1d_nan(a_1d_nan):
     return b
 
 
+@pytest.fixture
+def a_1d_with_zeros(a_with_zeros):
+    """Timeseries of a with zeros"""
+    return a_with_zeros.isel(lon=0, lat=0, drop=True)
+
+
 # weights
 @pytest.fixture
 def weights(a):
diff --git a/xskillscore/tests/test_metric_results_accurate.py b/xskillscore/tests/test_metric_results_accurate.py
index a00e3d4..54f1481 100644
--- a/xskillscore/tests/test_metric_results_accurate.py
+++ b/xskillscore/tests/test_metric_results_accurate.py
@@ -2,7 +2,12 @@ import numpy as np
 import pytest
 import sklearn.metrics
 from scipy.stats import pearsonr, spearmanr
-from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
+from sklearn.metrics import (
+    mean_absolute_error,
+    mean_absolute_percentage_error,
+    mean_squared_error,
+    r2_score,
+)
 
 import xskillscore as xs
 from xskillscore.core.deterministic import (
@@ -23,6 +28,11 @@ xs_skl_metrics = [
     (r2, r2_score),
     (mse, mean_squared_error),
     (mae, mean_absolute_error),
+    (mape, mean_absolute_percentage_error),
+]
+
+xs_skl_metrics_with_zeros = [
+    (mape, mean_absolute_percentage_error),
 ]
 
 xs_scipy_metrics = [
@@ -34,7 +44,6 @@ xs_scipy_metrics = [
 
 
 xs_np_metrics = [
-    (mape, lambda x, y: np.mean(np.abs((x - y) / x))),
     (me, lambda x, y: np.mean(x - y)),
     (smape, lambda x, y: 1 / len(x) * np.sum(np.abs(y - x) / (np.abs(x) + np.abs(y)))),
 ]
@@ -69,6 +78,15 @@ def test_xs_same_as_skl_same_name(a_1d, b_1d, request):
     assert np.allclose(actual, expected)
 
 
+@pytest.mark.parametrize("xs_skl_metrics", xs_skl_metrics_with_zeros)
+def test_xs_same_as_skl_with_zeros(a_1d_with_zeros, b_1d, xs_skl_metrics):
+    """Tests xskillscore metric is same as scikit-learn metric."""
+    xs_metric, skl_metric = xs_skl_metrics
+    actual = xs_metric(a_1d_with_zeros, b_1d, "time")
+    expected = skl_metric(a_1d_with_zeros, b_1d)
+    assert np.allclose(actual, expected)
+
+
 @pytest.mark.parametrize("xs_scipy_metrics", xs_scipy_metrics)
 def test_xs_same_as_scipy(a_1d, b_1d, xs_scipy_metrics):
     """Tests xskillscore metric is same as scipy metric."""