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,
>>> SimArray([1,2], "Mpc") + SimArray([1,2], "Msol"))
ValueError
>>> SimArray([1,2], "Mpc") + SimArray([1,2], "kpc")
SimArray([1.001, 1.002], "Mpc")
>>> 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
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')
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:
@SimArray.ufunc_rule(np.add)
@SimArray.ufunc_rule(np.subtract)
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
Defines a shallow wrapper around numpy.ndarray for extra functionality like unit-tracking.
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.
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.
x.__ge__(y) <==> x>=y
Same as self.transpose(), except that self is returned if self.ndim < 2.
Examples
>>> 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.
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:
@SimDict.getter
def X_copy(d) :
return d['X']
@SimDict.setter
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 (snap.dm, snap.star, 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.
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.
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 pynbody.family.Family 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. f.dm in place of the equivalent syntax f[pynbody.family.dm].
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.
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 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()).
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.
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 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.
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.
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.
Wraps the positions of the particles in the box to lie between [-boxsize/2, boxsize/2].
If no boxsize is specified, self.properties[“boxsize”] is used.
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.
Input
array_name - the name of the array to write
Optional Keywords
Returns True if the array or family array of given name is auto-derived (and therefore read-only).
Represent a sub-view of a SimSnap, initialized by specifying a slice. Arrays accessed through __getitem__ are automatically sub-viewed using the given slice.
Represents a subset of the simulation particles according to an index array.
Represents a one-family portion of a parent snap object
Loads a file using the appropriate class, returning a SimSnap instance.
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.
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 * units.cm**(-2,3)
In the last example, either a tuple describing a fraction or a Fraction instance (from the standard python module fractions) is acceptable.
To convert one unit to another, use the ratio member function:
>>> units.Msol.ratio(units.kg)
1.99e30
>>> (units.Msol / units.kpc**3).ratio(units.m_p/units.cm**3)
4.04e-8
If the units cannot be converted, a UnitsException is raised:
>>> units.Msol.ratio(units.kpc)
UnitsException
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)
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.
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")
1.000
>>> ((q**2)/(4*math.pi*epsilon0*units.m**2)).ratio("N")
2.31e-28
Base class for units. To instantiate a unit, call the pynbody.units.Unit() factory function.
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)
250.0
>>> Unit("1 Mpc").ratio("Msol")
UnitsException: not convertible
>>> Unit("1 Mpc").ratio("Msol", kg=25.0, m=50.0)
3.1028701506345152e-08
Class factory for units. Given a string s, creates a Unit object.
for example:
“1.e30 kg”
“kpc**2”
“26.2 m s**-1”
Returns a decorator to create a function which auto-converts input to given units.
Usage:
@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
Various utility routines used internally by pynbody.
Open a file, determining from the filename whether to use gzip decompression
Open a file for reading, returning also the (decompressed) file size
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.
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.
Return a slice s3 with the property that ar[s1][s2] == ar[s3]
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].
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.
Returns True if a1 and a2 are numpy views pointing to the exact same underlying data; False otherwise.
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.
Returns the index to the first element in array which satisfies array[index]>=find. The array must be sorted in ascending order.
Given an array ar, return nbins+1 monotonically increasing bin edges such that the number of items in each bin is approximately equal.
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
(http://tog.acm.org/resources/GraphicsGems/gemsiii/rand_rotation.c)
Incomplete gamma function accepting complex z, based on algorithm given in numerical recipes (3rd ed)
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.
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
‘sum’: sum over the outputs