Running rpy2 in parallel using multiprocessing raises weird exception that cannot be caught

Question

So this is a problem that I have not been able to solve, and neither do I know of a good way to make a MCVE out of. Essentially, it has been briefly discussed here, but as the comments show, there was some disagreement, and the final verdict is still out. Hence I am posting a similar question again, hoping to get a better answer.

Background

I have sensor data from a couple of thousand sensors, that I get every minute. My interest lies in forecasting the data. For this I am using the ARIMA family of forecasting models. Long story short, after discussion with the rest of my research group, we decided to use the Arima function available in the R package forecast, instead of the statsmodels implementation of the same.

Problem Definition

Since, I have data from a few thousand sensors, for which I would like to at least analyse a whole week's worth of data (to begin with), and since a week has 7 days, I have 7 times the number of sensors data with me. Essentially a some 14k sensor-day combinations. Finding the best ARIMA order (which minimizes BIC) and forecasting the next day of week data takes about 1 minute for each sensor-day combination. Which means upwards of 11 days to just process one week data on a single core!

This is obviously a waste, when I have 15 more cores just idling away the whole time. So, obviously, this is a problem for parallel processing. Note that each sensor-day combination does not influence any other sensor-day combination. Also, the rest of my code is fairly well profiled, and optimized.

Issue

The issue is that I get this weird error that I cannot catch anywhere. Here is the error reproduced:

Exception in thread Thread-3:
Traceback (most recent call last):
  File "/home/kartik/miniconda3/lib/python3.5/threading.py", line 914, in _bootstrap_inner
    self.run()
  File "/home/kartik/miniconda3/lib/python3.5/threading.py", line 862, in run
    self._target(*self._args, **self._kwargs)
  File "/home/kartik/miniconda3/lib/python3.5/multiprocessing/pool.py", line 429, in _handle_results
    task = get()
  File "/home/kartik/miniconda3/lib/python3.5/multiprocessing/connection.py", line 251, in recv
    return ForkingPickler.loads(buf.getbuffer())
  File "/home/kartik/miniconda3/lib/python3.5/site-packages/rpy2/robjects/robject.py", line 55, in _reduce_robjectmixin
    rinterface_level=rinterface_factory(rdumps, rtypeof)
ValueError: Mismatch between the serialized object and the expected R type (expected 6 but got 24)

Here are a few characteristics of this error that I have discovered:

It is raised in the rpy2 package
It has something to do with Thread 3. Since Python is zero indexed, I am guessing this is the fourth thread. Therefore, 4x6 = 24, which adds up to the numbers shown in the final error statement
rpy2 is being used in only one place in my code where it might have to recode returned values into Python types. Protecting that line in try: ... except: ... does not catch that exception
The exception is not raised when I ditch the multiprocessing and call the function within a loop
The exception does not crash the program, just suspends it forever (till I Ctrl+C it into terminating)
All that I tried till now, have had no effect in resolving the error

Things Tried

I have tried everything from extreme procedural coding, with functions to deal with the least cases (that is only one function to be called in parallel), to extreme encapsulation, where the executable block in the if __name__ == '__main__': calls a function which reads in the data, does the necessary grouping, then passes the groups to another function, which imports multiprocessing and calls another function in parallel, which imports the processing module that imports rpy2, and passes the data to the Arima function in R.

Basically, it doesn't matter if rpy2 is called and initialized deep inside function nests, such that it has no idea another instance might be initialized, or if it is called and initialized once, globally, the error is raised if multiprocessing is involved.

Pseudo Code

Here is an attempt to present at least some basic pseudo code such that the error might be reproduced.

import numpy as np
import pandas as pd

def arima_select(y, order):
    from rpy2 import robjects as ro
    from rpy2.robjects.packages import importr
    from rpy2.robjects import pandas2ri
    pandas2ri.activate()
    forecast = importr('forecast')

    res = forecast.Arima(y, order=ro.FloatVector(order))
    return res

def arima_wrapper(data):
    data = data[['tstamp', 'val']]
    data.set_index('tstamp', inplace=True)
    return arima_select(data, (1,1,1))

def applyParallel(groups, func):
    from multiprocessing import Pool, cpu_count
    with Pool(cpu_count()) as p:
        ret_list = p.map(func, [group for _, group in groups])
    return pd.concat(ret_list, keys=[name for name, _ in groups])

def wrapper():
    df = pd.read_csv('file.csv', parse_dates=[1], infer_datetime_format=True)
    df['day'] = df['tstamp'].dt.day
    res = applyParallel(df.groupby(['sensor', 'day']), arima_wrapper)
    print(res)

Obviously, the above code can be encapsulated further, but I think it should reproduce the error quite accurately.

Data Sample

Here is the output of print(data.head(6)) when placed immediately below data.set_index('tstamp', inplace=True) in arima_wrapper from the pseudo code above:

Or alternatively, data for a sensor, for a whole week can be generated simply with:

def data_gen(start_day):
    r = pd.Series(pd.date_range('2016-09-{}'.format(str(start_day)),
                                periods=24*60, freq='T'),
                  name='tstamp')
    d = pd.Series(np.random.randint(10, 80, 1440), name='val')
    s = pd.Series(['sensor1']*1440, name='sensor')
    return pd.concat([s, r, d], axis=1)
df = pd.concat([data_gen(day) for day in range(1,8)], ignore_index=True)

Observations and Questions

The first observation is that this error is only raised when multiprocessing is involved, not when the function (arima_wrapper) is called in a loop. Therefore, it must be associated somehow with multiprocessing issues. R is not very multiprocess friendly, but when written in the way shown in the pseudo code, each instance of R should not know about the existence of the other instances.

The way the pseudo code is structured, there must be an initialization of rpy2 for each call inside the multiple subprocesses spawned by multiprocessing. If that were true, each instance of rpy2 should have spawned its own instance of R, which should just execute one function, and terminate. That would not raise any errors, because it would be similar to the single threaded operation. Is my understanding here accurate, or am I completely or partially missing the point?

Were all instances of rpy2 to somehow share an instance of R, then I might reasonably expect the error. What is true: is R shared among all instances of rpy2, or is there an instance of R for each instance of rpy2?

How might this issue be overcome?

Since SO hates question threads with multiple questions in them, I will prioritize my questions such that partial answers will be accepted. Here is my priority list:

How might this issue be overcome? A working code example that does not raise the issue will be accepted as answer even if it does not answer any other question, provided no other answer does better, or was posted earlier.
Is my understanding of Python imports accurate, or am I missing the point about multiple instances of R? If I am wrong, how should I edit the import statements such that a new instance is created within each subprocess? Answers to this question are likely to point me towards a probable solution, and will be accepted, provided no answer does better, or was posted earlier
Is R shared among all instances of rpy2 or is there an instance of R for each instance of rpy2? Answers to this question will be accepted only if they lead to a resolution of the problem.

Kartik · Accepted Answer

The following changes to the arima_select function in the pesudo code presented in the question work:

import numpy as np
import pandas as pd
from rpy2 import rinterface as ri

ri.initr()

def arima_select(y, order):

    def rimport(packname):
        as_environment = ri.baseenv['as.environment']
        require = ri.baseenv['require']
        require(ri.StrSexpVector([packname]),
                quiet = ri.BoolSexpVector((True, )))
        packname = ri.StrSexpVector(['package:' + str(packname)])
        pack_env = as_environment(packname)
        return pack_env

    frcst = rimport("forecast")
    args = (('y', ri.FloatSexpVector(y)),
            ('order', ri.FloatSexpVector(order)),
            ('include.constant', ri.StrSexpVector(const)))
    return frcst['Arima'].rcall(args, ri.globalenv)

Keeping the rest of the pseudo code the same. Note that I have since optimized the code further, and it does not require all the functions presented in the question. Basically, the following is necessary and sufficient:

import numpy as np
import pandas as pd
from rpy2 import rinterface as ri

ri.initr()

def arima(y, order=(1,1,1)):
    # This is the same as arima_select above, just renamed to arima
    ...

def applyParallel(groups, func):
    from multiprocessing import Pool, cpu_count
    with Pool(cpu_count(), worker_init) as p:
        ret_list = p.map(func, [group for _, group in groups])
    return pd.concat(ret_list, keys=[name for name, _ in groups])

def main():
    # Create your df in your favorite way:
    def data_gen(start_day):
        r = pd.Series(pd.date_range('2016-09-{}'.format(str(start_day)),
                                    periods=24*60, freq='T'),
                      name='tstamp')
        d = pd.Series(np.random.randint(10, 80, 1440), name='val')
        s = pd.Series(['sensor1']*1440, name='sensor')
        return pd.concat([s, r, d], axis=1)
    df = pd.concat([data_gen(day) for day in range(1,8)], ignore_index=True)

    applyParallel(df.groupby(['sensor', pd.Grouper(key='tstamp', freq='D')]),
                  arima) # Note one may use partial from functools to pass order to arima

Note that I also do not call arima directly from applyParallel since my goal is to find the best model for the given series (for a sensor and day). I use a function arima_wrapper to iterate through the order combinations, and call arima at each iteration.