Essential Framework Modules


Defines a shallow wrapper around numpy.ndarray for extra functionality like unit-tracking.

For most purposes, the differences between numpy.ndarray and array.SimArray are not important. However, when units are specified (by setting the units attribute), the behaviour is slightly different. In particular,

  • it becomes impossible to add or subtract arrays with incompatible dimensions
>>> SimArray([1,2], "Mpc") + SimArray([1,2], "Msol"))
  • addition or subtraction causes auto-unit conversion. For example
>>> SimArray([1,2], "Mpc") + SimArray([1,2], "kpc")
SimArray([1.001, 1.002], "Mpc")
  • Note that in this context the left value takes precedence in specifying the return units, so that reversing the order of the operation here would return results in kpc.
  • If only one of the arrays specifies a Unit, no checking occurs and the unit of the returned array is assumed to be the same as the one specified input unit.
  • Powers to single integer or rational powers will maintain unit tracking. Powers to float or other powers will not be able to do so.
>>> SimArray([1,2],"Msol Mpc**-3")**2
SimArray([1, 4], 'Msol**2 Mpc**-6')
>>> SimArray([1,2],"Msol Mpc**-3")**(1,3)
SimArray([ 1.,1.26], 'Msol**1/3 Mpc**-1')

Syntax above mirrors syntax in units module, where a length-two tuple can represent a rational number, in this case one third.

>>> SimArray([1.,2], "Msol Mpc**-3")**0.333
SimArray([ 1.,1.26])  # Lost track of units

Getting the array in specified units

Given an array, you can convert it in-place into units of your own chosing:

>>> x = SimArray([1,2], "Msol")
>>> x.convert_units('kg')
>>> print x
SimArray([  1.99e+30,   3.98e+30], 'kg')

Or you can leave the original array alone and get a copy in different units, correctly converted:

>>> x = SimArray([1,2], "Msol")
>>> print x.in_units("kg")
SimArray([  1.99e+30,   3.98e+30], 'kg')
>>> print x
SimArray([1,2], "Msol")

If the SimArray was created by a SimSnap (which is most likely), it has a pointer into the SimSnap’s properties so that the cosmological context is automatically fetched. For example, comoving -> physical conversions are correctly achieved:

>>> f = pynbody.load("fname")
>>> f['pos']
SimArray([[ 0.05419805, -0.0646539 , -0.15700017],
         [ 0.05169899, -0.06193341, -0.14475258],
         [ 0.05672406, -0.06384531, -0.15909944],
         [ 0.0723075 , -0.07650762, -0.07657281],
         [ 0.07166634, -0.07453796, -0.08020873],
         [ 0.07165282, -0.07468577, -0.08020587]], '2.86e+04 kpc a')
>>> f['pos'].convert_units('kpc')
>>> f['pos']
SimArray([[ 1548.51403101, -1847.2525312 , -4485.71463308],
         [ 1477.1124212 , -1769.52421398, -4135.78377699],
         [ 1620.68592366, -1824.15000686, -4545.69387564],
         [ 2065.9264273 , -2185.92982874, -2187.79225915],
         [ 2047.60759667, -2129.6537339 , -2291.6758134 ],
         [ 2047.2214441 , -2133.87693163, -2291.59406997]], 'kpc')

Specifying rules for ufunc’s

In general, it’s not possible to infer what the output units from a given ufunc should be. While numpy built-in ufuncs should be handled OK, other ufuncs will need their output units defined (otherwise a numpy.ndarray will be returned instead of our custom type.)

To do this, decorate a function with SimArray.ufunc_rule(ufunc). The function you define should take the same number of parameters as the ufunc. These will be the input parameters of the ufunc. You should return the correct units for the output, or raise units.UnitsException (in the latter case, the return array will be made into a numpy.ndarray.)

For example, here is the code for the correct addition/subtraction handler:

def _consistent_units(a,b) :

    # This will be called whenever the standard numpy ufuncs np.add
    # or np.subtract are called with parameters a,b.

    # You should always be ready for the inputs to have no units.

    a_units = a.units if hasattr(a, 'units') else None
    b_units = b.units if hasattr(b, 'units') else None

    # Now do the logic. If we're adding incompatible units,
    # we want just to get a plain numpy array out. If we only
    # know the units of one of the arrays, we assume the output
    # is in those units.

    if a_units is not None and b_units is not None :
        if a_units==b_units :
            return a_units
        else :
            raise units.UnitsException("Incompatible units")

    elif a_units is not None :
        return a_units
    else :
        return b_units
class pynbody.array.SimArray[source]

Defines a shallow wrapper around numpy.ndarray for extra functionality like unit-tracking.


Provides the basemost SimArray that an IndexedSimArray is based on.


Set the units for this array by performing dimensional analysis on the supplied unit and referring to the units of the original file


Set the units for this array by performing dimensional analysis on the default dimensions for the array.


Retun a copy of this array expressed in the units specified in the parameter file.

in_units(new_unit, **context_overrides)[source]

Return a copy of this array expressed relative to an alternative unit.


Convert units of this array in-place. Note that if this is a sub-view, the entire base array will be converted.


Write this array to disk according to the standard method associated with its base file. This is equivalent to calling

>>> sim.gas.write_array('array')

in the case of writing out the array ‘array’ for the gas particle family. See the description of pynbody.snapshot.SimSnap.write_array() for options.

pynbody.array.wrapper_function(self, other, comparison_op=<slot wrapper '__ge__' of 'numpy.ndarray' objects>)[source]

x.__ge__(y) <==> x>=y


Same as self.transpose(), except that self is returned if self.ndim < 2.


>>> x = np.array([[1.,2.],[3.,4.]])
>>> x
array([[ 1.,  2.],
       [ 3.,  4.]])
>>> x.T
array([[ 1.,  3.],
       [ 2.,  4.]])
>>> x = np.array([1.,2.,3.,4.])
>>> x
array([ 1.,  2.,  3.,  4.])
>>> x.T
array([ 1.,  2.,  3.,  4.])

x.__ge__(y) <==> x>=y


This submodule defines an augmented dictionary class (SimDict) for SimSnap properties where entries need to be managed e.g. for defining default entries, or for ensuring consistency between equivalent properties like redshift and scalefactor.

By default, a SimDict automatically converts between redshift (‘z’) and scalefactor (‘a’) and implements default entries for cosmological values listed in the [default-cosmology] section of the pynbody configuration files.

Adding further properties

To add further properties use the SimDict.getter and SimDict.setter decorators. For instance, to add a property ‘X_copy’ which just reflects the value of the property ‘X’, you would use the following code:

def X_copy(d) :
    return d['X']

def X_copy(d, value) :
    d['X'] = value


This module defines the Family class which represents families of particles (e.g. dm, gas, star particles). New Family objects are automatically registered so that snapshots can use them in the normal syntax (,, etc).

In practice the easiest way to make use of the flexibility this module provides is through adding more families of particles in your config.ini.[source]

Returns a list of the names of all particle families. If with_aliases is True, include aliases in the list., create=False)[source]

Returns a family corresponding to the specified string. If the family does not exist and create is False, raises ValueError. If the family does not exist and create is True, an appropriate object is instantiated, registered and returned.


This module implements the SimSnap class which manages and stores snapshot data. It also implements the SubSnap class (and relatives) which represent different views of an existing SimSnap.

class pynbody.snapshot.SimSnap[source]

The class for managing simulation snapshots.

For most purposes, SimSnaps should be initialized through load() or new().

For a basic tutorial explaining how to load a file as a SimSnap see tutorials/data_access.

Getting arrays or subsnaps

Once a SimSnap object f is instantiated, it can be used in various ways. The most common operation is to access something with the code f[x]. Depending on the type of x, various behaviours result:

  • If x is a string, the array named by x is returned. If no such array exists, the framework attempts to load or derive an array of that name (in that order). If this is unsuccessful, a KeyError is raised.

  • If x is a python slice (e.g. f[5:100:3]) or an array of integers (e.g. f[[1,5,100,200]]) a subsnap containing only the mentioned particles is returned.

    See tutorials/data_access for more information.

  • If x is a numpy array of booleans, it is interpreted as a mask and a subsnap containing only those particles for which x[i] is True. This means that f[condition] is a shortcut for f[np.where(condition)].

  • If x is a pynbody.filt.Filter object, a subsnap containing only the particles which pass the filter condition is returned.

    See tutorials/data_access for more information.

  • If x is a object, a subsnap containing only the particles in that family is returned. In practice for most code it is more convenient to write e.g. in place of the equivalent syntax f[].

Getting metadata

The property filename gives the filename of a snapshot.

There is also a properties dictionary which contains further metadata about the snapshot. See Subsnaps.


Return the particle families which have representitives in this SimSnap. The families are ordered by their appearance in the snapshot.


Return the directly accessible array names (in memory)


Returns True if the array name is accessible (in memory)


Returns a list of the actual arrays in memory


Returns a list of tuples describing the array names and their contents in memory

get(key, alternative=None)[source]

Standard python get method, returns self[key] if key in self else alternative


Returns True if the array name is accessible (in memory) for at least one family


Returns a list of arrays which can be lazy-loaded from an auxiliary file.


Returns a list of arrays which can be lazy-evaluated.


Returns a list of all arrays that can be either lazy-evaluated or lazy loaded from an auxiliary file.


Return list of arrays which are not accessible from this view, but can be accessed from family-specific sub-views.

If fam is not None, only those keys applying to the specific family will be returned (equivalent to self.fam.keys()).


Returns true if other is a subview of self


Returns true if self is a subview of other


The original SimSnap from which this view is derived (potentially self)

get_index_list(relative_to, of_particles=None)[source]

Get a list specifying the index of the particles in this view relative to the ancestor relative_to, such that relative_to[get_index_list(relative_to)]==self.

intersect(other, op=<function intersect1d at 0x105032500>)[source]

Returns the set intersection of this simulation view with another view of the same simulation


Returns the set union of this simulation view with another view of the same simulation


Returns the set difference of this simulation view with another view of the same simulation


Return a dictionary containing a (scalefactor) and h (Hubble constant in canonical units) for this snapshot, ready for passing into unit conversion functions.

set_units_system(velocity=None, distance=None, mass=None, temperature=None)[source]

Set the unit system for the snapshot by specifying any or all of velocity, distance, mass and temperature units. The units can be given as strings or as pynbody Unit objects.

If any of the units are not specified and a previous file_units_system does not exist, the defaults are used.


Converts all arrays’units to be consistent with the units of the original file.

physical_units(distance='kpc', velocity='km s^-1', mass='Msol', persistent=True)[source]

Converts all array’s units to be consistent with the distance, velocity, mass basis units specified.

Base units can be specified using keywords.

Optional Keywords:

distance: string (default = ‘kpc’)

velocity: string (default = ‘km s^-1’)

mass: string (default = ‘Msol’)

persistent: boolean (default = True); apply units change to future lazy-loaded arrays if True


Given a unit (or string) dimensions, returns a unit with the same physical dimensions which is in the unit schema of the current file.

halos(*args, **kwargs)[source]

Tries to instantiate a halo catalogue object for the given snapshot, using the first available method (as defined in the configuration files).


Tries to construct a bridge function between this SimSnap and another one.

This function calls pynbody.bridge.bridge_factory(). For more information see bridge-tutorial, or the reference documentation for pynbody.bridge.


Rotates the snapshot about the current x-axis by ‘angle’ degrees.


Rotates the snapshot about the current y-axis by ‘angle’ degrees.


Rotates the snapshot about the current z-axis by ‘angle’ degrees.

wrap(boxsize=None, convention='center')[source]

Wraps the positions of the particles in the box to lie between [-boxsize/2, boxsize/2].

If no boxsize is specified,[“boxsize”] is used.

write_array(array_name, fam=None, overwrite=False, **kwargs)[source]

Write out the array with the specified name.

Some of the functionality is available via the pynbody.array.SimArray.write() method, which calls the present function with appropriate arguments.


array_name - the name of the array to write

Optional Keywords

fam (None) - Write out only one family; or provide a list to
write out a set of families.
is_derived_array(name, fam=None)[source]

Returns True if the array or family array of given name is auto-derived (and therefore read-only).

If the named array is auto-derived, this destroys the link so that the array becomes editable but no longer auto-updates.


Calculate the mean by mass of the specified array.

class pynbody.snapshot.SubSnap(base, _slice)[source]

Represent a sub-view of a SimSnap, initialized by specifying a slice. Arrays accessed through __getitem__ are automatically sub-viewed using the given slice.


Return the units on disk for a quantity with the specified dimensions

class pynbody.snapshot.IndexedSubSnap(base, index_array)[source]

Represents a subset of the simulation particles according to an index array.

class pynbody.snapshot.FamilySubSnap(base, fam)[source]

Represents a one-family portion of a parent snap object

pynbody.snapshot.load(filename, *args, **kwargs)[source]

Loads a file using the appropriate class, returning a SimSnap instance., order=None, **families)[source]

Create a blank SimSnap, with the specified number of particles.

Position, velocity and mass arrays are created and filled with zeros.

By default all particles are taken to be dark matter. To specify otherwise, pass in keyword arguments specifying the number of particles for each family, e.g.

f = new(dm=50, star=25, gas=25)

The order in which the different families appear in the snapshot is unspecified unless you add an ‘order’ argument:

f = new(dm=50, star=25, gas=25, order=’star,gas,dm’)

guarantees the stars, then gas, then dark matter particles appear in sqeuence.


The pynbody units module consists of a set of classes for tracking units.

It relates closely to the array module, which defines an extension to numpy arrays which carries unit information.

Making units

Units are generated and used at various points through the pynbody framework. Quite often the functions where users interact with units simply accept strings.

You can also make units yourself in two ways. Either you can create a string, and instantiate a Unit like this:

>>> units.Unit("Msol kpc**-3")
>>> units.Unit("2.1e12 m_p cm**-2/3")

Or you can do it within python, using the named Unit objects

>>> units.Msol * units.kpc**-3
>>> 2.1e12 * units.m_p ***(-2,3)

In the last example, either a tuple describing a fraction or a Fraction instance (from the standard python module fractions) is acceptable.

Getting conversion ratios

To convert one unit to another, use the ratio member function:

>>> units.Msol.ratio(
>>> (units.Msol / units.kpc**3).ratio(units.m_p/**3)

If the units cannot be converted, a UnitsException is raised:

>>> units.Msol.ratio(units.kpc)

Specifying numerical values

Sometimes it’s necessary to specify a numerical value in the course of a conversion. For instance, consider a comoving distance; this can be specified in pynbody units as follows:

>>> comoving_kpc = units.kpc * units.a

where units.a represents the scalefactor. We can attempt to convert this to a physical distance as follows

>>> comoving_kpc.ratio(units.kpc)

but this fails, throwing a UnitsException. On the other hand we can specify a value for the scalefactor when we request the conversion

>>> comoving_kpc.ratio(units.kpc, a=0.5)

and the conversion completes with the expected result. The units module also defines units.h for the dimensionless hubble constant, which can be used similarly. By default, all conversions happening within a specific simulation context should pass in values for a and h as a matter of routine.

Any IrreducibleUnit (see below) can have a value specified in this way, but a and h are envisaged to be the most useful applications.

Defining new base units

The units module is fully extensible: you can define and name your own units which then integrate with all the standard functions.

litre = units.NamedUnit("litre",0.001*units.m**3)
gallon = units.NamedUnit("gallon",0.004546*units.m**3)
gallon.ratio(litre) # 4.546
(units.pc**3).ratio(litre) # 2.94e52

You can even define completely new dimensions.

V = units.IrreducibleUnit("V") # define a volt
C = units.NamedUnit("C", units.J/V) # define a coulomb
q = units.NamedUnit("q", 1.60217646e-19*C) # elementary charge
F = units.NamedUnit("F", C/V) # Farad
epsilon0 = 8.85418e-12 *F/units.m
>>> (q*V).ratio("eV")
>>> ((q**2)/(4*math.pi*epsilon0*units.m**2)).ratio("N")
class pynbody.units.UnitBase[source]

Base class for units. To instantiate a unit, call the pynbody.units.Unit() factory function.

ratio(other, **substitutions)[source]

Get the conversion ratio between this Unit and another specified unit.

Keyword arguments, if specified, give numerical substitutions for the named unit. This is most useful for specifying values for cosmological quantities like ‘a’ and ‘h’, but can also be used for any IrreducibleUnit.

>>> Unit("1 Mpc a").ratio("kpc", a=0.25)
>>> Unit("1 Mpc").ratio("Msol")
UnitsException: not convertible
>>> Unit("1 Mpc").ratio("Msol", kg=25.0, m=50.0)
in_units(*a, **kw)[source]

Alias for ratio


Return a unit equivalent to this one (may be identical) but expressed in terms of the currently defined IrreducibleUnit instances.


Class factory for units. Given a string s, creates a Unit object.

The string format is:
[<scale>] [<unit_name>][**<rational_power>] [[<unit_name>] ... ]

for example:

“1.e30 kg”


“26.2 m s**-1”

pynbody.units.takes_arg_in_units(*args, **orig_kwargs)[source]

Returns a decorator to create a function which auto-converts input to given units.


@takes_arg_in_units((2, "Msol"), (1, "kpc"), ("blob", "erg"))
def my_function(arg0, arg1, arg2, blob=22) :
   print "Arg 2 is",arg2,"Msol"
   print "Arg 1 is",arg1,"kpc"
   print "blob is",blob,"ergs"
>>> My_function(22, "1.e30 kg", 23, blob="12 J")
Input 3 is 0.5 Msol
Input 2 is 23 kpc

Returns True if the specified object has a meaningful units attribute


Returns True if the specified object has a meaningful units attribute


Returns True if the specified object represents a unit


Returns True if the specified object is itself a unit or otherwise exposes unit information


Various utility routines used internally by pynbody.

pynbody.util.open_(filename, *args)[source]

Open a file, determining from the filename whether to use gzip decompression

pynbody.util.open_with_size(filename, *args)[source]

Open a file for reading, returning also the (decompressed) file size

pynbody.util.intersect_slices(s1, s2, array_length=None)[source]

Given two python slices s1 and s2, return a new slice which will extract the data of an array d which is in both d[s1] and d[s2].

Note that it may not be possible to do this without information on the length of the array referred to, hence all slices with end-relative indexes are first converted into begin-relative indexes. This means that the slice returned may be specific to the length specified.

pynbody.util.relative_slice(s_relative_to, s)[source]

Given a slice s, return a slice s_prime with the property that array[s_relative_to][s_prime] == array[s]. Clearly this will not be possible for arbitrarily chosen s_relative_to and s, but it should be possible for s=intersect_slices(s_relative_to, s_any) which is the use case envisioned here (and used by SubSim). This code currently does not work with end-relative (i.e. negative) start or stop positions.

pynbody.util.chained_slice(s1, s2)[source]

Return a slice s3 with the property that ar[s1][s2] == ar[s3]

pynbody.util.index_before_slice(s, index)[source]

Return an index array new_index with the property that, for a slice s (start, stop and step all positive), ar[s][index] == ar[new_index].

pynbody.util.concatenate_indexing(i1, i2)[source]

Given either a numpy array or slice for both i1 and i2, return either a numpy array or slice i3 with the property that

ar[i3] == ar[i1][i2].

As a convenience, if i2 is None, i1 is returned


Given either an array or slice, return len(ar[sl_or_ar]) for any array ar which is large enough that the slice does not overrun it.

pynbody.util.arrays_are_same(a1, a2)[source]

Returns True if a1 and a2 are numpy views pointing to the exact same underlying data; False otherwise.

pynbody.util.set_array_if_not_same(a_store, a_in, index=None)[source]

This routine checks whether a_store and a_in ultimately point to the same buffer; if not, the contents of a_in are copied into a_store.

pynbody.util.index_of_first(array, find)[source]

Returns the index to the first element in array which satisfies array[index]>=find. The array must be sorted in ascending order.

pynbody.util.equipartition(ar, nbins, min=None, max=None)[source]

Given an array ar, return nbins+1 monotonically increasing bin edges such that the number of items in each bin is approximately equal.

pynbody.util.bisect(left, right, f, epsilon=None, eta=0, verbose=False, niter_max=200)[source]

Finds the value x such that f(x)=0 for a monotonically increasing function f, using a binary search.

The search stops when either the bounding domain is smaller than epsilon (by default 10^-7 times the original region) OR a value f(x) is found such that |f(x)|<eta (by default eta=0, so this criterion is never satisfied).


A simple Gauss-Jordan matrix inverter. This is provided so that matrices of fractions can be inverted (numpy linalg converts everything to floats first.)

Don’t use on large matrices – it’s slow!

Based on public domain code by Jarno Elonen.


A simple replacement for numpy linalg matrix inverse which handles fractions exactly. Not suitable for large matrices!


Return a random rotation matrix (Haar measure for 3x3 case), using fast algorithm from Graphics Gems III



Strip the .gz ending off a string

pynbody.util.gamma_inc(a, z, eps=3e-07)[source]

Incomplete gamma function accepting complex z, based on algorithm given in numerical recipes (3rd ed)

pynbody.util.threadsafe_inline(*args, **kwargs)[source]

When scipy.weave.inline is called, it may trigger a compile. We only want one compilation to be going on at once, otherwise nasty race conditions arise. This function wraps scipy.weave.inline to be thread-safe.

pynbody.util.parallel(p_args=[0], threads=4, reduce='interleave')[source]

Return a function decorator which makes a function execute in parallel.

p_args: a list of integers specifying which arguments will be array-like. These will be sliced up and processed over the number of threads specified.

reduce: specifies how to reduce the output, and can take one of three values:

‘none’: return None

‘interleave’: return an array of the size of the input arrays,
with the elements returned to their correct place

‘sum’: sum over the outputs