Repository: IAMconsortium/pyam

Issue:
`interpolate` does not follow other return conventions
We recently had a private communication regarding confusion with the `interpolate` function. The user was trying to do something like:

```
df_interp = df.filter().interpolate()
```

The returned object was *not* an interpolated dataframe. We advised the user as follows:

```
alldata = pyam.IamDataFrame(data='iamc15_scenario_data_world_r1.1.xlsx')
df = alldata.filter(variable='Emissions|CO2',scenario='SSP1-19') # save the filtered dataframe explicitly
df.interpolate(2063)  # interpolate the filtered dataframe
```

I would argue this breaks our implicit datamodel of dataframe methods returning a copy unless `inplace=True` is provided. 

@danielhuppmann can you provide your thoughts?

While making such a change would break existing API, I think there is a strong argument for all functions to conform to our (from `pandas`) data model.

Discussion:
gidden: seems like test failures are related to iiasa queries!
danielhuppmann: > seems like test failures are related to iiasa queries!

already pinged @peterkolp and @fonfon about this...
danielhuppmann: @gidden, is this related to #371 - not quite sure what the bug to be resolved is?
znicholls: Looks like @danielhuppmann has the review under control. My only thought: if you cast to wide format first, could you get an extra performance boost by using panda’s or numpy’s inbuilt Interpolation rather than looping over all the times?
coroa: > My only thought: if you cast to wide format first, could you get an extra performance boost by using panda’s or numpy’s inbuilt Interpolation rather than looping over all the times?

Agreed.

```python
if not islistable(time):
    time = [time]
pyam.IamDataFrame(
    idf.timeseries().reindex(columns=np.sort(idf.year + time)).interpolate(method='slinear', axis=1),
    meta=idf.meta
)
```
has proven to beat pyam's implementation by heaps and bounds.

danielhuppmann: Indeed, it would be great to switch to a more performant implementation of the interpolation!

One issue to keep in mind, though, is that [pandas.DataFrame.interpolate](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.interpolate.html) fills in all missing values in the dataframe, not just the values for one new column. That's why I didn't use it at the time. And [numpy.interp](https://numpy.org/doc/stable/reference/generated/numpy.interp.html) assumes that the values are constant outside the timeseries domain, which is also not a valid general assumption.
znicholls: > is that [pandas.DataFrame.interpolate](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.interpolate.html) fills in all missing values in the dataframe, not just the values for one new column

Yes you'd have to choose the subset you want to interpolate first if you didn't just want to interpolate everything. 

Another thought: you could just also change the API so the user just specifies the new time points they desire (including any which are already present) rather than just adding **extra** time points alone.

So, say you had data for 2010, 2020 and 2030 and you wanted to interpolate to have 2010, 2015, 2020, 2025 and 2030. 

In the current implementation you'd do

```python
df.interpolate([2015, 2025])
```

But we could just make it so it's

```python
df.interpolate([2010, 2015, 2020, 2025, 2030])
```

And just leave all the calculations to pandas. That interface would mean that interpolating would always leave you with a nan-free dataset.

> And [numpy.interp](https://numpy.org/doc/stable/reference/generated/numpy.interp.html) assumes that the values are constant outside the timeseries domain

[scipy interp1d](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d) will extrapolate if you want with the `fill_value` argument.

p.s. @coroa 

> heaps and bounds

was that a typo?
coroa: > p.s. @coroa
> 
> > heaps and bounds
> 
> was that a typo?



> # Please confirm that this PR has done the following:
> 
>     * [x]  Tests Added
> 
>     * [x]  Documentation Added
> 
>     * [ ]  Description in RELEASE_NOTES.md Added
> 
> 
> ## Adding to RELEASE_NOTES.md (remove section after adding to RELEASE_NOTES.md)
> 
> Please add a single line in the release notes similar to the following:
> 
> ```
> - (#XX)[http://link-to-pr.com] Added feature which does something
> ```
> 
> # Description of PR
> 
> This PR fixes a bug when there are no values to interpolate (all scenarios already have this value) and allows users to past a list-like input for interpolation of multiple time points.



> > is that [pandas.DataFrame.interpolate](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.interpolate.html) fills in all missing values in the dataframe, not just the values for one new column
> 
> Yes you'd have to choose the subset you want to interpolate first if you didn't just want to interpolate everything.
> 
> Another thought: you could just also change the API so the user just specifies the new time points they desire (including any which are already present) rather than just adding **extra** time points alone.
> 
> So, say you had data for 2010, 2020 and 2030 and you wanted to interpolate to have 2010, 2015, 2020, 2025 and 2030.
> 
> In the current implementation you'd do
> 
> ```python
> df.interpolate([2015, 2025])
> ```
> 
> But we could just make it so it's
> 
> ```python
> df.interpolate([2010, 2015, 2020, 2025, 2030])
> ```
> 
> And just leave all the calculations to pandas. That interface would mean that interpolating would always leave you with a nan-free dataset.
> 
> > And [numpy.interp](https://numpy.org/doc/stable/reference/generated/numpy.interp.html) assumes that the values are constant outside the timeseries domain
> 
> [scipy interp1d](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d) will extrapolate if you want with the `fill_value` argument.
>

`method='slinear'` of pandas is indeed passed to scipy's `interp1d` (together with all extra `**kwargs`) and its default is to raise a `ValueError`, unless `fill_value` is set to f.ex. `'extrapolate'.

> p.s. @coroa
> 
> > heaps and bounds
> 
> was that a typo?

something like that :)


coroa: Otherwise, something along the lines of (untested)
```python
df = self.timeseries()
new_values = df.apply(lambda s: pd.Series(scipy.interpolate.interp1d(s.index, s.values, method='slinear', **kwargs)(time), time), axis=1)
for year in time:
    if year in df:
        df[year].fillna(new_values[year], inplace=True)
    else:
        i = bisect(df.columns, year)
        df.insert(i, year, new_values[year])
return pyam.IamDataFrame(df, meta=self.meta)
```
should work fine.
gidden: Thanks for the discussion all! I'll take a look at the suggestions, update, and flag for final review.
gidden: Ok folks, I have updated implementation as attached in this notebook. Let me know what you think - it should maintain prior behavior (only updating columns asked for). 
[interpolate_implementation (1).pdf](https://github.com/IAMconsortium/pyam/files/5540789/interpolate_implementation.1.pdf)

gidden: Ok... this took me a bit further down a rabbit hole than I expected. The basic approach follows what @coroa suggested earlier, and includes:

1. pivot to wide-format once
2. interpolate using `pandas.DataFrame.interpolate`, passing along kwargs
3. include back only time values that were asked for (keeping old behavior)
4. return a copy by default (in the future)

There were a few issues I ran into that I will raise to others. Namely an issue for @danielhuppmann / @Rlamboll regarding interpolation with `extra_col`s and an issue with `datetime.datetime` for @znicholls. I will ping on the specific lines.

@coroa happy if you take a look at the final version here too.
gidden: Ok @danielhuppmann - should be good for final review now. I must admit, I am a bit surprised at the rather minor speedups (@coroa can you confirm these numbers match your observations?). But in general, this PR was meant to support additional input, not optimization, so no need to block it there.
danielhuppmann: resolved merge conflicts

Patch:
diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md
index 83bb415..ea75eca 100644
--- a/RELEASE_NOTES.md
+++ b/RELEASE_NOTES.md
@@ -5,6 +5,7 @@
 - [#461](https://github.com/IAMconsortium/pyam/pull/461) Add list of authors to repo and docs pages
 - [#459](https://github.com/IAMconsortium/pyam/pull/459) Add a `get_variable_components()` function to retrieve or join variable components
 - [#458](https://github.com/IAMconsortium/pyam/pull/458) Enable `Path` for IamDataFrame initialization 
+- [#456](https://github.com/IAMconsortium/pyam/pull/456) Speed up interpolation and support interpolation for multiple time-points
 - [#454](https://github.com/IAMconsortium/pyam/pull/454) Enable dimensionless units and fix `info()` if IamDataFrame is empty
 - [#451](https://github.com/IAMconsortium/pyam/pull/451) Fix unit conversions from C to CO2eq
 - [#450](https://github.com/IAMconsortium/pyam/pull/450) Defer logging set-up to when the first logging message is generated
diff --git a/pyam/core.py b/pyam/core.py
index 5fb4b26..942b21a 100755
--- a/pyam/core.py
+++ b/pyam/core.py
@@ -507,34 +507,71 @@ class IamDataFrame(object):
         df = df.unstack(level=columns, fill_value=fill_value)
         return df
 
-    def interpolate(self, time):
-        """Interpolate missing values in timeseries (linear interpolation)
+    def interpolate(self, time, inplace=None, **kwargs):
+        """Interpolate missing values in the timeseries data
+
+        This method uses :meth:`pandas.DataFrame.interpolate`,
+        which applies linear interpolation by default
 
         Parameters
         ----------
-        time : int or datetime
+        time : int or datetime, or list-like thereof
              Year or :class:`datetime.datetime` to be interpolated.
              This must match the datetime/year format of `self`.
+        inplace : bool, optional
+            if True, do operation inplace and return None
+        kwargs
+            passed to :meth:`pandas.DataFrame.interpolate`
         """
-        if self.time_col == 'year' and not isinstance(time, int):
+        ##
+        # TODO remove, and add kwarg inplace=False in release >= 0.10
+        if inplace is None:
+            deprecation_warning(
+                'Behavior of `interpolate` will change to `inplace=False` '
+                'as default in a future release. Set the kwarg explicitly '
+                'to avoid this warning. Use `inplace=True` to keep current '
+                'behavior.'
+            )
+            inplace = True
+        ##
+
+        # setup
+        ret = self.copy() if not inplace else self
+        interp_kwargs = dict(method='slinear', axis=1)
+        interp_kwargs.update(kwargs)
+        time = list(time) if islistable(time) else [time]
+        # TODO - have to explicitly cast to numpy datetime to sort later,
+        # could enforce as we do for year below
+        if self.time_col == 'time':
+            time = list(map(np.datetime64, time))
+        elif not all(isinstance(x, int) for x in time):
             raise ValueError(
-                'The `time` argument `{}` is not an integer'.format(time)
+                'The `time` argument `{}` contains non-integers'.format(time)
             )
+        old_cols = list(ret[ret.time_col].unique())
+        columns = np.sort(np.unique(old_cols + time))
+
+        # calculate a separate dataframe with full interpolation
+        df = ret.timeseries()
+        newdf = (
+            df
+            .reindex(columns=columns)
+            .interpolate(**interp_kwargs)
+        )
 
-        # get data in wide format, drop rows where value is already defined
-        df = self.timeseries()
-        if time in df.columns:
-            df = df[np.isnan(df[time])]
+        # replace only columns asked for
+        for col in time:
+            df[col] = newdf[col]
 
-        # apply fill_series, re-add time dimension to index, set series name
-        _values = df.apply(fill_series, raw=False, axis=1, time=time).dropna()
-        _values.index = append_index_level(
-            index=_values.index, codes=[0] * len(_values),
-            level=[time], name=self.time_col, order=self._data.index.names)
-        _values.name = 'value'
+        # replace underlying data object
+        # TODO naming time_col could be done in timeseries()
+        df.columns.name = ret.time_col
+        df = df.stack()  # long-data to pd.Series
+        df.name = 'value'
+        ret._data = df.sort_index()
 
-        # append interpolated values to `_data` and sort index
-        self._data = self._data.append(_values).sort_index()
+        if not inplace:
+            return ret
 
     def swap_time_for_year(self, inplace=False):
         """Convert the `time` column to `year`.


Test patch:
diff --git a/tests/test_core.py b/tests/test_core.py
index ca48099..b7825d0 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -631,6 +631,8 @@ def test_interpolate(test_pd_df):
     df.interpolate(2007)
     obs = df.filter(year=2007).data['value'].reset_index(drop=True)
     exp = pd.Series([3, 1.5, 4], name='value')
+    print(obs)
+    print(exp)
     pd.testing.assert_series_equal(obs, exp)
 
     # redo the interpolation and check that no duplicates are added
@@ -641,10 +643,53 @@ def test_interpolate(test_pd_df):
     assert all([True if isstr(i) else ~np.isnan(i) for i in df.data.foo])
 
 
+def test_interpolate_time_exists(test_df_year):
+    df = test_df_year
+    df.interpolate(2005)
+    obs = df.filter(year=2005).data['value'].reset_index(drop=True)
+    exp = pd.Series([1.0, 0.5, 2.0], name='value')
+    pd.testing.assert_series_equal(obs, exp)
+
+
+def test_interpolate_with_list(test_df_year):
+    df = test_df_year
+    df.interpolate([2007, 2008])
+    obs = df.filter(year=[2007, 2008]).data['value'].reset_index(drop=True)
+    exp = pd.Series([3, 4, 1.5, 2, 4, 5], name='value')
+    pd.testing.assert_series_equal(obs, exp)
+
+
+def test_interpolate_full_example():
+    cols = ['model_a', 'scen_a', 'World']
+    df = IamDataFrame(pd.DataFrame([
+        cols + ['all', 'EJ/yr', 0, 1, 6., 10],
+        cols + ['last', 'EJ/yr', 0, 0.5, 3, np.nan],
+        cols + ['first', 'EJ/yr', 0, np.nan, 2, 7],
+        cols + ['middle', 'EJ/yr', 0, 1, np.nan, 7],
+        cols + ['first two', 'EJ/yr', 0, np.nan, np.nan, 7],
+        cols + ['last two', 'EJ/yr', 0, 1, np.nan, np.nan],
+    ], columns=IAMC_IDX + [2000, 2005, 2010, 2017],
+    ))
+    exp = IamDataFrame(pd.DataFrame([
+        cols + ['all', 'EJ/yr', 0, 1, 6., 7.142857, 10],
+        cols + ['last', 'EJ/yr', 0, 0.5, 3, np.nan, np.nan],
+        cols + ['first', 'EJ/yr', 0, 1., 2, 3.428571, 7],
+        cols + ['middle', 'EJ/yr', 0, 1, np.nan, 4.5, 7],
+        cols + ['first two', 'EJ/yr', 0, 2.058824, np.nan, 4.941176, 7],
+        cols + ['last two', 'EJ/yr', 0, 1, np.nan, np.nan, np.nan],
+    ], columns=IAMC_IDX + [2000, 2005, 2010, 2012, 2017],
+    ))
+    obs = df.interpolate([2005, 2012], inplace=False)
+    assert_iamframe_equal(obs, exp)
+
+
 def test_interpolate_extra_cols():
-    # check hat interpolation with non-matching extra_cols has no effect (#351)
+    # check that interpolation with non-matching extra_cols has no effect
+    # (#351)
     EXTRA_COL_DF = pd.DataFrame([
         ['foo', 2005, 1],
+        ['foo', 2010, 2],
+        ['bar', 2005, 2],
         ['bar', 2010, 3],
     ],
         columns=['extra_col', 'year', 'value'],
@@ -656,8 +701,11 @@ def test_interpolate_extra_cols():
     df2 = df.copy()
     df2.interpolate(2007)
 
-    # assert that interpolation didn't change any data
-    assert_iamframe_equal(df, df2)
+    # interpolate should work as if extra_cols is in the _data index
+    assert_iamframe_equal(df, df2.filter(year=2007, keep=False))
+    obs = df2.filter(year=2007)['value']
+    exp = pd.Series([2.4, 1.4], name='value')
+    pd.testing.assert_series_equal(obs, exp)
 
 
 def test_interpolate_datetimes(test_df):