[1]:

import numpy as np
import scipp as sc

Tips, tricks, and anti-patterns#

Choose dimensions wisely#

A good choice of dimensions for representing data goes a long way in making working with Scipp efficient. Consider, e.g., data gathered from detector pixels at certain time intervals. We could represent it as

[2]:

npix = 100
ntime = 10
data = sc.zeros(dims=['pixel', 'time'], shape=[npix, ntime])
data

[2]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (8.06 KB)
    • (pixel: 100, time: 10)
      float64
      𝟙
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

For irregularly spaced detectors this may well be the correct or only choice. If however the pixels are actually forming a regular 2-D image sensor we should probably prefer

[3]:

nx = 10
ny = npix // nx
data = sc.zeros(dims=['y', 'x', 'time'], shape=[ny, nx, ntime])
data

[3]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (8.06 KB)
    • (y: 10, x: 10, time: 10)
      float64
      𝟙
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([[[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]])

With this layout we can naturally perform slices, access neighboring pixel rows or columns, or sum over rows or columns.

Choose dimension order wisely#

In principle the order of dimensions in Scipp can be arbitrary since operations transpose automatically based on dimension labels. In practice however a bad choice of dimension order can lead to performance bottlenecks. This is most obvious when slicing multi-dimensional variables or arrays, where slicing any but the outer dimension yields a slice with gaps between data values, i.e., a very inefficient memory layout. If an application requires slicing (directly or indirectly, e.g., in groupby operations) predominantly for a certain dimension, this dimension should be made the outermost dimension. For example, for a stack of images the best choice would typically be

[4]:

nimage = 13
images = sc.zeros(
    dims=['image', 'y', 'x'],
    shape=[
        nimage,
        ny,
        nx,
    ],
)
images

[4]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (10.41 KB)
    • (image: 13, y: 10, x: 10)
      float64
      𝟙
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([[[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       ...,

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]])

Slices such as

[5]:

images['image', 3]

[5]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (1.03 KB out of 10.41 KB)
    • (y: 10, x: 10)
      float64
      𝟙
      0.0, 0.0, ..., 0.0, 0.0
      Values:
      array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

will then have data for all pixels in a contiguous chunk of memory. Note that in Scipp the first listed dimension in dims is always the outermost dimension (NumPy’s default).

Avoid loops#

With Scipp, just like with NumPy or Matlab, loops such as for-loops should be avoided. Loops typically lead to many small slices or many small array objects and rapidly lead to very inefficient code. If we encounter the need for a loop in a workflow using Scipp we should try and take a step back to understand how it can be avoided. Some tips to do this include:

Use slicing with “shifts”#

When access to neighbor slices is required, replace

[6]:

for i in range(len(images.values) - 1):
    images['image', i] -= images['image', i + 1]

with

[7]:

images['image', :-1] -= images['image', 1:]

Note that at this point NumPy provides more powerful functions such as numpy.roll. Scipp’s toolset for such purposes is not fully developed yet.

Seek advice from NumPy#

There is a huge amount of information available for NumPy, e.g., on stackoverflow. We can profit in two ways from this. In some cases, the same techniques can be applied to Scipp variables or data arrays, since mechanisms such as slicing and basic operations are very similar. In other cases, e.g., when functionality is not available in Scipp yet, we can resort to processing the raw array accessible through the values property:

[8]:

var = sc.arange('x', 10.0)
var.values = np.roll(var.values, 2)
var

[8]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (336 Bytes)
    • (x: 10)
      float64
      𝟙
      8.0, 9.0, ..., 6.0, 7.0
      Values:
      array([8., 9., 0., 1., 2., 3., 4., 5., 6., 7.])

The values property can also be used as the out argument that many NumPy functions support:

[9]:

np.exp(var.values, out=var.values)
var

[9]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (336 Bytes)
    • (x: 10)
      float64
      𝟙
      2980.958, 8103.084, ..., 403.429, 1096.633
      Values:
      array([2.98095799e+03, 8.10308393e+03, 1.00000000e+00, 2.71828183e+00,
       7.38905610e+00, 2.00855369e+01, 5.45981500e+01, 1.48413159e+02,
       4.03428793e+02, 1.09663316e+03])

WARNING

When applying NumPy functions to the values directly we lose handling of units and variances, so this should be used with care.

Use helper dimensions or reshaped data#

Some operations may be difficult to implement without a loop in a certain data layout. If this layout cannot be changed globally, we can still change it temporarily for a certain operation. Even if this requires a copy it may still be faster and more concise than implementing the operation with a loop. For example, we can sum neighboring elements by temporarily reshaping with a helper dimension using fold:

[10]:

var = sc.arange('x', 12.0)
var.fold('x', sizes={'x': 4, 'neighbors': 3}).sum('neighbors')

[10]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (288 Bytes)
    • (x: 4)
      float64
      𝟙
      3.0, 12.0, 21.0, 30.0
      Values:
      array([ 3., 12., 21., 30.])

Note that fold returns a view, i.e., the operation is performed without making a copy of the underlying data buffers. The companion operation of fold is flatten, which provides the reverse operation (see the section below for more details).

Use in-place operations#

Allocating memory or copying data is an expensive process and may even be the dominant factor for overall application performance, apart from loading large amounts of data from disk. Therefore, it pays off the avoid copies where possible.

Scipp provides two mechanisms for this, in-place arithmetic operators such as +=, and out-arguments similar to what NumPy provides. Examples:

[11]:

var = var * 2.0  # makes a copy
var *= 2.0  # in-place (faster)

[12]:

var = sc.sqrt(var)  # makes a copy
var = sc.sqrt(var, out=var)  # in-place (faster)

Note that in-place operations cannot be used if a broadcast is required or a dtype change happens, since in-place operations may only change the data contained in a variable.

Reshaping data#

The shape of a Variable or a DataArray can be modified using the fold and flatten functions. Below are a few examples to illustrate how they work.

Folding#

In a nutshell, the fold operation increases the number of dimensions of the data. We begin with a two-dimensional variable:

[13]:

N = 4
M = 3
var = sc.array(dims=['x', 'y'], values=np.random.random([N, M]))
sc.show(var)
dims=('x', 'y'), shape=(4, 3), unit=dimensionless, variances=Falsevalues xy

We then fold the x dimension into two new dimensions a and b:

[14]:

folded_var = sc.fold(var, dim='x', sizes={'a': 2, 'b': 2})
sc.show(folded_var)
dims=('a', 'b', 'y'), shape=(2, 2, 3), unit=dimensionless, variances=Falsevalues a by

The result is a three-dimensional variable with dimensions (a, b, y).

A DataArray with coordinates can also be folded:

[15]:

x = sc.array(dims=['x'], values=np.arange(N))
y = sc.array(dims=['y'], values=np.arange(M))
da = sc.DataArray(data=var, coords={'x': x, 'y': y})
sc.show(da)
(dims=('x', 'y'), shape=(4, 3), unit=dimensionless, variances=False)values xy xx(dims=('x',), shape=(4,), unit=dimensionless, variances=False)values x yy(dims=('y',), shape=(3,), unit=dimensionless, variances=False)values y
[16]:

folded_da = sc.fold(da, dim='x', sizes={'a': 2, 'b': 2})
sc.show(folded_da)
(dims=('a', 'b', 'y'), shape=(2, 2, 3), unit=dimensionless, variances=False)values a by xx(dims=('a', 'b'), shape=(2, 2), unit=dimensionless, variances=False)values a b yy(dims=('y',), shape=(3,), unit=dimensionless, variances=False)values y

Note that the dimensions of the x coordinate have changed from (x,) to (a, b), but the coordinate name has not changed.

Flattening#

The inverse of the fold operation is flatten. This is analogous to NumPy’s flatten method. By default, all dimensions of the input are flattened to a single dimension, whose name is provided by the to argument:

[17]:

flat_da = sc.flatten(da, to='z')
sc.show(flat_da)
(dims=('z',), shape=(12,), unit=dimensionless, variances=False)values z xx(dims=('z',), shape=(12,), unit=dimensionless, variances=False)values z yy(dims=('z',), shape=(12,), unit=dimensionless, variances=False)values z

It is however possible to only flatten selected dimensions, using the dims argument. For example, we can flatten the a and b dimensions of our previously folded (three-dimensional) data to recover a two-dimensional array.

[18]:

flat_ab = sc.flatten(folded_da, dims=['a', 'b'], to='time')
sc.show(flat_ab)
(dims=('time', 'y'), shape=(4, 3), unit=dimensionless, variances=False)values timey xx(dims=('time',), shape=(4,), unit=dimensionless, variances=False)values time yy(dims=('y',), shape=(3,), unit=dimensionless, variances=False)values y

Stacking: concatenating into a new dimension#

Another very common operation is combining multiple arrays together. In NumPy, stack is used to combine arrays along a new dimension, while concatenate is used to combine arrays along an existing dimension.

Because of its labeled dimensions, Scipp can achieve both operations using concat.

For example, giving concat a dim which is not found in the inputs will stack the arrays along a new dimension:

[19]:

stacked = sc.concat([da, da], dim='z')
sc.show(stacked)
(dims=('z', 'x', 'y'), shape=(2, 4, 3), unit=dimensionless, variances=False)values z xy xx(dims=('x',), shape=(4,), unit=dimensionless, variances=False)values x yy(dims=('y',), shape=(3,), unit=dimensionless, variances=False)values y

Improving performance: Allocators and HugePages#

Overview#

Scipp frequently needs to allocate or deallocate large chunks of memory. In many cases this is a major contribution to the (lack of) performance of Scipp (or any other) applications. Aside from minimizing the number of allocations (which is not always easy or possible), there are two simple ways that have each shown to independently yield up to 5% or 10% performance improvements for certain workloads.

Note

The following sections are for Linux only. Similar options may be available on other platforms, please consult the appropriate documentation.

Using the TBB allocation proxy#

Scipp uses TBB for parallelization. TBB also provides allocators, but Scipp is not using them by default. It is possible to automatically do so by preloading libtbbmalloc_proxy.so. For more context refer to the oneTBB documentation.

Note

Make sure to install Scipp via conda since Scipp’s PyPI package does not ship the required libtbbmalloc_proxy.so.

Python#

Use

$ LD_PRELOAD=libtbbmalloc_proxy.so.2 python

to launch Python or

$ LD_PRELOAD=libtbbmalloc_proxy.so.2 python myapp.py

to run a Python script. Note that the .2 suffix may be different depending on the actual version of TBB.

Jupyter#

Environment variables can be configured for a specific Jupyter kernel. First, obtain a list of kernels:

$ jupyter kernelspec list

The list may include, e.g., /home/<username>/mambaforge/envs/myenv/share/jupyter/kernels/python3. Consider making a copy of this folder, e.g., to /home/<username>/mambaforge/envs/myenv/share/jupyter/kernels/python3-tbbmalloc_proxy. Next, edit the kernel.json file in this folder, updating the “display_name” and adding an “env” section:

{
 "argv": [
  "/home/<username>/mambaforge/envs/myenv/bin/python",
  "-m",
  "ipykernel_launcher",
  "-f",
  "{connection_file}"
 ],
 "display_name": "Python 3 (ipykernel) with tbbmalloc_proxy",
 "language": "python",
 "env":{
   "LD_PRELOAD":"libtbbmalloc_proxy.so.2"
 }
}

Notes:

  • The .2 suffix may be different depending on the actual version of TBB.

  • File content aside from the “env” section may be different depending on your kernel. The above is an example and not a suggestion to change other settings.

Using HugePages#

HugePages are a Linux feature that may decrease the overhead of page handling and translating virtual addresses when dealing with large amounts of memory. This can be beneficial for Scipp since it frequently allocates and accesses a lot of memory.

Warning

HugePages is a complex topic and this does not aim to be a comprehensive guide. In particular, details on your system may be different from what is described here, please carefully check any commands given below before executing them, as they require root privileges. Please contact your system administrator for more information.

Using explicit HugePages#

HugePages can be reserved explicitly, in contrast to the transparent HugePages described below. First, check the number of configured HugePages:

$ cat /sys/kernel/mm/hugepages/hugepages-2048kB/nr_hugepages

By default this is 0 on many systems. We can change this at system runtime (as root),

# echo 1024 > /sys/kernel/mm/hugepages/hugepages-2048kB/nr_hugepages

where 1024 is the numfer of 2 MByte pages to allocate. This will allocate 2 GByte of memory as HugePages. You will likely need more than this — adjust this number depending on your system’s memory and the memory requirements of your application.

Unlike transparent HugePages, explicit HugePages are not used by default by an application. The TBB allocator can be configured to use HugePages by setting the TBB_MALLOC_USE_HUGE_PAGES environment variable. For example, we can extend the “env” section of the Jupyter kernel configuration above to

{
 "argv": [
  "/home/<username>/mambaforge/envs/myenv/bin/python",
  "-m",
  "ipykernel_launcher",
  "-f",
  "{connection_file}"
 ],
 "display_name": "Python 3 (ipykernel) with tbbmalloc_proxy",
 "language": "python",
 "env":{
   "LD_PRELOAD":"libtbbmalloc_proxy.so.2",
   "TBB_MALLOC_USE_HUGE_PAGES":"1"
 }
}

Using HugePages transparently#

The Linux kernel can be configured to transparently use HugePages for all allocations. We can check the setting using

$ cat /sys/kernel/mm/transparent_hugepage/enabled
always [madvise] never

The madvise setting, which is the default on some systems, is not sufficient for Scipp to use HugePages, even when using libtbbmalloc_proxy.so. Instead, we need to set the transparent_hugepage kernel parameter to always:

# echo always > /sys/kernel/mm/transparent_hugepage/enabled

Notes:

  • This enables hugepages globally and may adversely affect other applications.

  • This does not require use of libtbbmalloc_proxy.so, but the two can be used together.

Checking if it works#

Run

$ grep HugePages /proc/meminfo

to check the current use or

$ watch grep HugePages /proc/meminfo

to monitor the use while running your application. When using explicit HugePages the size in HugePages_Total should correspond to the setting made above. Now run your application:

  • When using transparent HugePages the size in AnonHugePages should increase. Note that NumPy appears to be using transparent HugePages as well (even when support is set to madvise instead of always), make sure you are not misled by this.

  • When using explicit HugePages the size in HugePages_Free should decrease.