title: Why does scipy.interpolate.griddata hang when used with multiprocessing and OpenBLAS?
body: <p>I have random samples of 3D data (x,y,z) where x and y are spatial dimensions and z is a function of its position on the plane: z = f(x,y). I want to evaluate this function over a regular grid using interpolated values from the samples.</p>
  
  <p>Each set of samples is in a different file and needs to be interpolated separately (each file is a different point in time). In order to speed up the processing of these files, I want to use the multiprocessing (mp) module to do them in parallel using a mp pool. </p>
  
  <p>However, when I try to use mp.griddata inside a function called by a mp pool, the process hangs when it gets to the griddata function. No problems when the same function is executed serially. Also this only happens with certain numpy/scipy builds. When my numpy.show_config() is the following: </p>
  
  <pre><code>blas_info:
      libraries = ['blas']
      library_dirs = ['/usr/lib']
      language = f77
  lapack_info:
      libraries = ['lapack']
      library_dirs = ['/usr/lib']
      language = f77
  atlas_threads_info:
    NOT AVAILABLE
  blas_opt_info:
      libraries = ['blas']
      library_dirs = ['/usr/lib']
      language = f77
      define_macros = [('NO_ATLAS_INFO', 1)]
  atlas_blas_threads_info:
    NOT AVAILABLE
  lapack_opt_info:
      libraries = ['lapack', 'blas']
      library_dirs = ['/usr/lib']
      language = f77
      define_macros = [('NO_ATLAS_INFO', 1)]
  atlas_info:
    NOT AVAILABLE
  lapack_mkl_info:
    NOT AVAILABLE
  blas_mkl_info:
    NOT AVAILABLE
  atlas_blas_info:
    NOT AVAILABLE
  mkl_info:
    NOT AVAILABLE
  </code></pre>
  
  <p>everything works fine.</p>
  
  <p>but when the config is:</p>
  
  <pre><code>lapack_info:
      libraries = ['openblas']
      library_dirs = ['/u/ccor/local/lib']
      language = f77
  atlas_threads_info:
      libraries = ['openblas']
      library_dirs = ['/u/ccor/local/lib']
      define_macros = [('ATLAS_WITHOUT_LAPACK', None)]
      language = c
      include_dirs = ['/usr/include']
  blas_opt_info:
      libraries = ['openblas']
      library_dirs = ['/u/ccor/local/lib']
      define_macros = [('ATLAS_INFO', '"\\"None\\""')]
      language = c
      include_dirs = ['/usr/include']
  atlas_blas_threads_info:
      libraries = ['openblas']
      library_dirs = ['/u/ccor/local/lib']
      define_macros = [('ATLAS_INFO', '"\\"None\\""')]
      language = c
      include_dirs = ['/usr/include']
  lapack_opt_info:
      libraries = ['openblas', 'openblas']
      library_dirs = ['/u/ccor/local/lib']
      define_macros = [('ATLAS_WITHOUT_LAPACK', None)]
      language = f77
      include_dirs = ['/usr/include']
  lapack_mkl_info:
    NOT AVAILABLE
  blas_mkl_info:
    NOT AVAILABLE
  mkl_info:
    NOT AVAILABLE
  </code></pre>
  
  <p>the child processes hang. </p>
  
  <p>The second configuration was created by trying to optimize my BLAS libraries for a machine with many cores by compiling LibOpenBLAS and numpy by hand. In both cases I am in a virtualenv using python 2.7.3.</p>
  
  <p>Also, the hanging only happens when griddata is used with 'linear' or 'cubic' interpolation; 'nearest' works fine. </p>
  
  <p>Here is some code trying to illustrate and isolate the problem. The data is randomly generated instead of reading from a file and the results are just accumulated arbitrarily as an example. </p>
  
  <pre><code>import multiprocessing as mp
  import numpy
  from scipy import interpolate
  
  def interpolate_grid(it, n_points, grid_dim):
  
      print 'interpolating: ', it
  
      numpy.random.seed(it)
  
      points = 2*numpy.random.random((n_points,2))-0.5
      values = numpy.random.random(n_points)
      x,y = numpy.mgrid[0:1:grid_dim*1j,0:1:grid_dim*1j]
      grid_pnts = numpy.vstack([x.flatten(), y.flatten()]).T
      interp_vals = interpolate.griddata(points, values, grid_pnts, 'cubic')
      print 'done interpolating: ', it
      return interp_vals
  
  def interpolation_test(n_runs = 10, n_points = 1000, grid_dim = 10, parallel = True):
  
      class Sum(object):
          def __init__(self, grid_dim):
              self.acc = numpy.zeros(grid_dim*grid_dim)
  
          def accumulate(self, val):
              self.acc += val
  
      acc = Sum(grid_dim)
  
      if parallel:
          pool = mp.Pool(mp.cpu_count())
          for i in xrange(n_runs):
              pool.apply_async(interpolate_grid, 
                              args = [i, n_points, grid_dim], 
                              callback = acc.accumulate)
          pool.close()
          pool.join()
  
      else:
          for i in xrange(n_runs):
              acc.acc += interpolate_grid(i, n_points, grid_dim)
  
      return acc.acc
  
  def main():
  
      ser_result = interpolation_test(parallel = False)
      par_result = interpolation_test(parallel = True)
  
      assert (numpy.abs(par_result-ser_result) &lt; 1e-14).all()
      print 'serial and parallel results the same'
  
  if __name__ == '__main__':
      main()
  </code></pre>
  
  <p>Edit: I found this on the <a href="http://scikit-learn.org/stable/install.html" rel="nofollow">scikit-learn installation page</a>: "Using OpenBLAS can give speedups in some scikit-learn modules, but it doesn’t play nicely with joblib/multiprocessing, so using it is not recommended unless you know what you’re doing."</p>
  
  <p>I found a workaround by using ATLAS instead of OpenBLAS:</p>
  
  <pre><code>aptitude install libatlas-base-dev
  ATLAS=/usr/lib
  pip install numpy
  </code></pre>
  
  <p>It's not quite as fast as OpenBLAS, but works with multiprocessing, which is more important for my application, and is much faster than the basic numpy installation. Does anyone who "knows what they are doing" care to comment on why OpenBLAS and multiprocessing "don't play well together"? Or how one might get them to play nice? The warning on the scikit-learn site seems to imply that it is possible.</p>