Index to pyEXP classes

pyEXP

Provides a collection of EXP tools for processing and analyzing simulation data using BFE techniques and MSSA.

How to get started

The main documentation is the many docstrings embedded in the code and a set of examples provided with the EXP source. We hope to provide a online reference guide in the future.

We recommend beginning with the example Python scripts and IPython notebooks and adapting them to your own needs. You can explore the available classes and member functions using the usual Python ``help’’ function. The classes are organized into seven submodules that are described briefly below. Run ‘help(pyEXP.xxxx) for each of the submodules below for more detailed usage info…

The available submodules

read: Read particle snapshots of various types. Currently EXP, Gadget, Tipsy, and Bonzai types are supported.
basis: Create and apply specific biorthogonal bases to generate coefficients from particle data and evaluate potential, density, and force fields
coefs: Classes for reading, passing, writing, converting, and querying coefficient sets
field: Create two- and three-dimension rectangular grids of fields for visualization
mssa: Tools to apply Multivariate Singular Spectrum Analysis (MSSA) to the coefficients computed using the ‘basis’ classes
edmd: Tools to apply the discrete Koopman operator analysis using the extended Dynamical Mode Decomposition (EDMD) algorith to approximate the Koopman operator. The use of the Koopman class echos that for mSSA. As of this point, this is NOT a recommended toolset. If you have success with this, please post a message.
util: Miscellaneous tools that support the others. Currently this include centering algorithms. While EXP has native methods for doing this, others will need to supply an estimated center

Example workflow

To provide some context, suppose you want to read some snapshots, make some coefficients, and then analyze them with MSSA. The class construction would go something like this:

Create a reader instance for your simulation, call it ‘reader’.

Create a basis designed to represent the a particular particle type. Star particles, perhaps, so let’s call it ‘disk’.

We then pass ‘reader’ to the createCoefficients member of ‘disk’ to get coefficients for your snapshots, called ‘coefs’

We might then want to explore dynamical patterns in these coefficients by passing ‘coefs’ to ‘expMSSA’. ‘expMSSA’ will return principal signals as an updated coefficient object, that we call ‘newcoefs’

‘newcoefs’ and ‘disk’ can then be passed to the FieldGenerator to provide density, potential, force fields, etc. for the each principal signal

This is only one example of many possible uses. There are many variants to this work flow, of course, and I expect that you will invent some interesting ones

The source code has some sample Python scripts and notebooks for a quick start (check the pyEXP directory).

Version information

Version: EXP 7.9.0 Repository URL: https://github.com/EXP-code/EXP GIT branch: SLboundaries GIT commit: 9b78e241ce712fd41c4b2430702ce9d6f2626078 Compile time: 2025-07-26 02:40:47 UTC

History and provenance

These EXP interface classes and the Python interface were written to fill a demand by collaborators for Python access to EXP. EXP itself represents 20 years of experience applying BFE techniques to both simulation and perturbation theory. Recent innovations include the application of MSSA to discover and characterize some key dynamical features in simulations that are hard to find `by eye’.

Please send comments, suggestions, and particularly good cookies to: mdw@umass.edu (Martin Weinberg)

ParticleReader class bindings

This collection of classes reads and converts your phase-space snapshots to iterable objects for generating basis coefficients. The particle fields themselves are not available in Python for inspection currently.

The available particle readers are:

PSPout The monolithic EXP phase-space snapshot format
PSPspl Like PSPout, but split into multiple file chunks

3. GadgetNative The original Gadget native format 4 GadgetHDF5 The newer HDF5 Gadget format 5. TipsyNative The original Tipsy format 6. TipsyXDR The original XDR Tipsy format 7. Bonsai This is the Bonsai varient of Tipsy files

We have a helper function, getReaders, to get a list to help you remember. Try: pyEXP.read.ParticleReader.getReaders()

Each reader can manage snapshots split into many files by parallel, per process writing. The reader classes takes a list of lists on input. The outer list represents individual snapshots and the inner list is all files belonging to a single snapshot. For unsplit snapshots, the inner list is a single file. We provide two helper functions that will make the list of lists from the lexical sort and pattern match on a single list of filenames of the form ‘prefix_xxxxx-yyyyy, where xxxxx is the integer index of the snapshot and yyyyy is the process index for a single snapshot. The value for ‘prefix’ is arbitrary. The delimiter between xxxxx and yyyyy may be user specified. See routines ‘parseFileList()’ and ‘parseStringList()’.

Once the ParticleReader instance is created, you can select the type of particle you want to read; each reader has different mnemonics for these, as you know. You can get the list of types using the GetTypes() member and select the desired type with the SelectType() member. You can also get the current time and number of particles for the selected type using the CurrentTime() and CurrentNumber() members, respectively. The main function of of this ParticleReader object in pyEXP is creating the coefficient expansion using Basis.createCoefficients. See help(pyEXP.basis)

BasisFactory class bindings

This module provides a factory class that will create biorthogonal bases from input YAML configuration files. Each basis can then be used to compute coefficients, provide field quantities such as forces and, together with the FieldGenerator, surfaces and fields for visualization.

Eight bases are currently implemented:

SphericalSL, the Sturm-Liouiville spherical basis;

Bessel, the classic spherical biorthogonal constructed from the eigenfunctions of the spherical Laplacian;

Cylindrical, created created by computing empirical orthogonal functions over a densely sampled SphericalSL basis;

FlatDisk, an EOF rotation of the finite Bessel basis;

CBDisk, the Clutton-Brock disk basis for testing;

Slab, a biorthogonal basis for a slab geometry with a finite finite vertical extent. The basis is constructed from direct solution of the Sturm-Liouville equation.

Cube, a periodic cube basis whose functions are the Cartesian eigenfunctions of the Cartesian Laplacian: sines and cosines.

FieldBasis, for computing user-provided quantities from a phase-space snapshot.

VelocityBasis, for computing the mean field velocity fields from a phase-space snapshot. This is a specialized version of FieldBasis.

Each of these bases take a YAML configuration file as input. These parameter lists are as subset of and have the same structure as thosed used by EXP. The factory and the individual constructors will check the parameters keys and warn of mismatches for safety. See the EXP documentation and the pyEXP examples for more detail. The first four bases are the most often used bi- orthogonal basis types used for computing the potential and forces from density distributions. Other biorthgonal bases in EXP but not in pyEXP include those for cubic and slab geometries and other special-purpose bases such as the Hernquist, Clutton-Brock sphere and two-dimensional disk basis. These will be made available in a future release if there is demand. Note that the Hernquist and Clutton-Brock spheres can be constructed using SphericalSL with a Hernquist of modified Plummer model as input. The FieldBasis and VelocityBasis are designed for producing summary data for post-production analysis (using mSSA or eDMD, for example) and for simulation cross-comparison.

The primary functions of these basis classes are:

To compute BFE coefficients from phase-space snapshots using the ParticleReader class. See help(pyEXP.read).

To evaluate the fields from the basis and a coefficient object. See help(pyEXP.coefs) and help(pyEXP.field).

To provide compact summary field data for post-production analysis. See help(pyEXP.basis.FieldBasis) and help(pyEXP.basis.VelocityBasis).

Introspection

The first two bases have a ‘cacheInfo(str)’ member that reports the parameters used to create the cached basis. This may be used to grab the parameters for creating a basis. Cache use ensures that your analyses are computed with the same bases used in a simulation or with the same basis used on previous pyEXP invocations. At this point, you must create the YAML configuration for the basis even if the basis is cached. This is a safety and consistency feature that may be relaxed in a future version.

Coefficient creation

The Basis class creates coefficients from phase space with two methods: ‘createFromReader()’ and ‘createFromArray()’. The first uses a ParticleReader, see help(pyEXP.read), and the second uses arrays of mass and 3d position vectors. Both methods take an optional center vector (default: 0, 0, 0). You may also register and an optional boolean functor used to select which particles to using the ‘setSelector(functor)’ member. An example functor would be defined in Python as follows:

def myFunctor(m, pos, vel, index):
ret = False # Default return value # some caculation with scalar mass, pos array, vel array and # integer index that sets ret to True if desired … return ret

If you are using ‘createFromArray()’, you will only have access to the mass and position vector. You may clear and turn off the selector using the ‘clrSelector()’ member.

The FieldBasis class requires a user-specified phase-space field functor that produces an list of quantities derived from the phase space for each particle. For example, to get a total velocity field, we could use:

def totalVelocity(m, pos, vel):
# Some caculation with scalar mass, pos array, vel array. # Total velocity for this example… return [(vel[0]**2 + vel[1]**2 + vel[2]**2)**0.5]

This function is registered with the FieldBasis using:

basis->addPSFunction(totalVelocity, [‘total velocity’])

The VelocityBasis is a FieldBasis that automatically sets the phase-space field functor to cylindrical or spherical velocities based on the ‘dof’ parameter. More on ‘dof’ below.

Scalablility

createFromArray() is a convenience method allows you to transform coordinates and preprocess phase space using your own methods and readers. Inside this method are three member functions calls that separately initialize, accumulate the coefficient contributions from the provided vectors, and finally construct and return the new coeffi- cient instance (Coefs). For scalability, we provide access to each of these three methods so that the phase space may be partitioned into any number of smaller pieces. These three members are: initFromArray(), addFromArray(), makeFromArray(). The initFromArray() is called once to begin the creation and the makeFromArray() method is called once to build the final set of coefficients. The addFromArray() may be called any number of times in between. For example, the addFromArray() call can be inside of a loop that iterates over any partition of phase space from your own pipeline. The underlying computation is identical to createFromArray(). However, access to the three underlying steps allows you to scale your phase-space processing to snapshots of any size. For reference, the createFromReader() method uses a producer-consumer pattern internally to provide scalability. These three methods allow you to provide the same pattern in your own pipeline.

Coordinate systems

Each basis is assigned a natural coordinate system for field evaluation as follows:

SphericalSL uses spherical coordinates

Cylindrical uses cylindrical coordinates

FlatDisk uses cylindrical coordinates

CBDisk uses cylindrical coordinates

Slab uses Cartesian coordinates

Cube uses Cartesian coordinates

FieldBasis and VelocityBasis provides two natural geometries for field evaluation: a two-dimensional (dof=2) polar disk and a three-dimensional (dof=3) spherical geometry that are chosen using the ‘dof’ parameter. These use cylindrical and spherical coordinates, respectively, by default.

These default choices may be overridden by passing a string argument to the ‘setFieldType()’ member. The argument is case insensitive and only distinguishing characters are necessary. E.g. for ‘Cylindrical’, the argument ‘cyl’ or even ‘cy’ is sufficient. The argument ‘c’ is clearly not enough.

Orbit integration

The IntegrateOrbits routine uses a fixed time step leap frog integrator to advance orbits from tinit to tfinal with time step h. The initial positions and velocities are supplied in an nx6 NumPy array. Tuples of the basis (a Basis instance) and coefficient database (a Coefs instance) for each component is supplied to IntegrateOrbtis as a list. Finally, the type of acceleration is an instance of the AccelFunc class. The acceleration at each time step is computed by setting a coefficient set in Basis and evaluating and accumulating the acceleration for each phase-space point. The coefficient are handled by implementing the evalcoefs() method of AccelFunc. We supply two implemented derived classes, AllTimeFunc and SingleTimeFunc. The first interpolates on the Coefs data base and installs the interpolated coefficients for the current time in the basis instance. The SingleTimeFunc interpolates on the Coefs data base for a single fixed time and sets the interpolated coefficients once at the beginning of the integration. This implementes a fixed potential model. AccelFunc can be inherited by a native Python class and the evalcoefs() may be implemented in Python and passed to IntegrateOrbits in the same way as a native C++ class.

Non-inertial frames of reference

Each component of a multiple component simulation may have its own expansion center. Orbit integration in the frame of reference of the expansion is accomplished by defining a moving frame of reference using the setNonInertial() call with either an array of n times and center positions (as an nx3 array) or by initializing with an EXP orient file.

We provide a member function, setNonInertialAccel(t), to estimate the frame acceleration at a given time. This may be useful for user-defined acceleration routines. This is automatically called default C++ evalcoefs() routine.

pyEXP.basis.IntegrateOrbits(tinit: float, tfinal: float, h: float, ps: numpy.ndarray[numpy.float64[m, n]], basiscoef: list[tuple[pyEXP.basis.Basis, pyEXP.coefs.Coefs]], func: pyEXP.basis.AccelFunc, nout: int = 0) → tuple[numpy.ndarray[numpy.float64[m, 1]], numpy.ndarray[numpy.float32]]

Compute particle orbits in gravitational field from the bases

Integrate a list of initial conditions from ‘tinit’ to ‘tfinal’ with a step size of ‘h’ using the list of basis and coefficient pairs. The step size will be adjusted to provide uniform sampling. Every step will be returned unless you provide an explicit value for ‘nout’, the number of desired output steps. In this case, the code will choose new set step size equal or smaller to the supplied step size with a stride to provide exactly ‘nout’ output times.

Parameters

tinitfloat: the intial time
tfinalfloat: the final time
hfloat: the integration step size
psnumpy.ndarray: an n x 6 table of phase-space initial conditions
bfelist(BasisCoef): a list of BFE coefficients used to generate the gravitational field
funcAccelFunctor: the force function
noutint: the number of output points, if specified

Returns

tuple(numpy.array, numpy.ndarray): time and phase-space arrays

Coefficient class bindings

These classes store, write, and provide an interface to coefficients and table data for use by the other pyEXP classes.

CoefStruct

The CoefStruct class is low-level structure that stores the data and metadata specific to each geometry. There are three groups of CoefStruct derived classes for biorthogonal basis coefficients, field data coeffients, and auxiliary table data. The biorthogonal classes are spherical (SphStruct), cylindrical (CylStruct), slab (SlabStruct), and cube (CubeStruct). The field classes cylindrical (CylFldStruct), and spherical (SphFldStruct). The table data is stored in TblStruct.

Instances of these structures represent individual times points and are created, maintained, and interfaced by the Coefs class. Access to the underlying data is provided to Python in case you need to change or rewrite the data for some reason. We have also provided a assign() member so that you can instaniate and load a coefficient structure using Python. To do this, use the constructor to make a blank instance, assign the dimensions and use assign() to create a data matrix with the supplied matrix or array. The dimen- sions are:

(lmax, nmax) for SphStruct

(mmax, nmax) for a CylStruct

(nmaxx, nmaxy, nmaxz) for a SlabStruct

(nmaxx, nmaxy, nmaxz) for a CubeStruct

(nfld, lmax, nmax) for a SphFldStruct

(nfld, mmax, nmax) for a CylFldStruct

(cols) for a TblStruct.

Coefs

The base class, ‘Coefs’, provides a factory reader that will create one of the derived coefficient classes, SphCoefs, CylCoefs, SlabCoefs, CubeCoefs, TblCoefs, SphFldCoefs, and CylFldCoefs, deducing the type from the input file. The input files may be EXP native or HDF5 cofficient files. Only biorthgonal basis coefficients have a native EXP type. The Basis factory, Basis::createCoefficients, will create set of coefficients from phase-space snapshots. See help(pyEXP.basis). Files which are not recognized as EXP coefficient files are assumed to be data files and are parsed by the TblCoefs class. The first column in data tables is interpreted as time and each successive column is interpreted as a new data field.

Once created, you may get a list of times, get the total gravitational power from biothogonal basis coefficients and general power from the field coefficients, and write a new HDF5 file. Their primary use is as a container object for MSSA (using expMSSA) and field visualization using the FieldGenerator class.

Updates

The expMSSA class will update the contribution to the coefficients specified by key from each eigen component to the reconstructed series. Unspecified coefficients series will not be updated and their original data will be intact. For visualization, the series data in a Coefs object may be zeroed using the ‘zerodata()’ member function prior to an expMSSA update. This allows one to include reconstructions that only include particular eigen components for the coefficients specified by key. Then, one can visualize only the updated fields using ‘FieldGenerator’. See help(pyEXP.mssa) and help(pyEXP.field) for more details.

Dataset indexing

Coefficients and other auxilliary data from simulations are stored and retrieved by their time field. Internally, these are floating fixed-point values truncated to 8 signficant figures so that they may be used as dictionary/map keys. The values in the time list, returned with the Times() member function contain the truncated values for reference.

Object lifetime

As in native Python, the memory for created objects persists until it is no longer referenced. For example, replacing the variable with a new set of coefficients will allow the memory to be deallocated if no other class instance holds a reference. Because coefficient sets can be large, the creation of a Coefs instance by the ‘factory’ will be passed to any other class that needs it. The MSSA class, expMSSA, will hold a reference to the the Coefs object passed on creation, and it update the values of the coefficients on reconstruction, without copying. If you want to keep the initial set without change, we have provided a ‘deepcopy()’ member that provides a byte-by-byte copy.

class pyEXP.coefs.CubeCoefs

Container for cube coefficients

PowerDim(self: pyEXP.coefs.CubeCoefs, d: str, min: int = 0, max: int = 2147483647) → numpy.ndarray[numpy.float64[m, n]]

Get power for the coefficient DB as a function of harmonic index for a given dimension. This Power() member is equivalent to PowerDim(‘x’).

Parameters

dchar: dimension for power summary; one of ‘x’, ‘y’, or ‘z’
minint: minimum index along requested dimension (default=0)
maxint: maximum index along requested dimension (default=max int)

Returns

numpy.ndarray: 2-dimensional numpy array containing the power

getAllCoefs(self: pyEXP.coefs.CubeCoefs) → numpy.ndarray[numpy.complex128]

Provide a 4-dimensional ndarray indexed by nx, ny, nz, and time indices.

Returns

numpy.ndarray: 4-dimensional numpy array containing the cube coefficients

setTensor(self: pyEXP.coefs.CubeCoefs, time: float, tensor: Eigen::Tensor<std::complex<double>, 3, 0, long>) → None

Enter and/or rewrite the coefficient tensor at the provided time

Parameters

timefloat: snapshot time corresponding to the the coefficient tensor mat : numpy.ndarray

the new coefficient array.

Returns

None

class pyEXP.coefs.CylCoefs

Container for cylindrical coefficients

EvenOddPower(self: pyEXP.coefs.CylCoefs, nodd: int = -1, min: int = 0, max: int = 2147483647) → tuple[numpy.ndarray[numpy.float64[m, n]], numpy.ndarray[numpy.float64[m, n]]]

Get cylindrical coefficient power separated into vertically even and odd contributions.

Parameters

noddint, default=-1: number of odd vertical modes to compute
minint, default=0: minimum time index for power calculation
maxint, default=inf: maximum time index for power calculation

Returns

numpy.ndarray: 2-dimensional numpy array containing the even and odd power coefficients

Notes

The default parameters (nodd<0) will query the YAML config for the value of ncylodd, but this can be provided as an argument if it is not explicitly set in your EXP::Cylinder configuration. If in doubt, use the default.

getAllCoefs(self: pyEXP.coefs.CylCoefs) → numpy.ndarray[numpy.complex128]

Provide a 3-dimensional ndarray indexed by azimuthal index, radial index, and time index

Returns

numpy.ndarray: 3-dimensional numpy array containing the cylindrical coefficients

setMatrix(self: pyEXP.coefs.CylCoefs, time: float, mat: numpy.ndarray[numpy.complex128[m, n]]) → None

Enter and/or rewrite the coefficient matrix at the provided time

Parameters

timefloat: snapshot time corresponding to the the coefficient matrix mat : numpy.ndarray

the new coefficient array.

Returns

None

class pyEXP.coefs.CylFldCoefs

Container for cylindrical field coefficients

getAllCoefs(self: pyEXP.coefs.CylFldCoefs) → numpy.ndarray[numpy.complex128]

Provide a 4-dimensional ndarray indexed by channel index, spherical index, radial index, and time index

Returns

numpy.ndarray: 4-dimensional numpy array containing the cylindrical coefficients

setMatrix(self: pyEXP.coefs.CylFldCoefs, time: float, mat: numpy.ndarray[numpy.complex128]) → None

Enter and/or rewrite the coefficient tensor at the provided time

Parameters

timefloat: snapshot time corresponding to the the coefficient matrix
matnumpy.ndarray: the new coefficient array.

Returns

None

class pyEXP.coefs.SlabCoefs

Container for cube coefficients

PowerDim(self: pyEXP.coefs.SlabCoefs, d: str, min: int = 0, max: int = 2147483647) → numpy.ndarray[numpy.float64[m, n]]

Get power for the coefficient DB as a function of harmonic index for a given dimension. This Power() member is equivalent to PowerDim(‘x’).

Parameters

dchar: dimension for power summary; one of ‘x’, ‘y’, or ‘z’
minint: minimum index along requested dimension (default=0)
maxint: maximum index along requested dimension (default=max int)

Returns

numpy.ndarray: 2-dimensional numpy array containing the power

getAllCoefs(self: pyEXP.coefs.SlabCoefs) → numpy.ndarray[numpy.complex128]

Provide a 4-dimensional ndarray indexed by nx, ny, nz, and time indices.

Returns

numpy.ndarray: 4-dimensional numpy array containing the slab coefficients

setTensor(self: pyEXP.coefs.SlabCoefs, time: float, tensor: Eigen::Tensor<std::complex<double>, 3, 0, long>) → None

Enter and/or rewrite the coefficient tensor at the provided time

Parameters

timefloat: snapshot time corresponding to the the coefficient tensor mat : numpy.ndarray

the new coefficient array.

Returns

None

class pyEXP.coefs.SphCoefs

Container for spherical coefficients

getAllCoefs(self: pyEXP.coefs.SphCoefs) → numpy.ndarray[numpy.complex128]

Provide a 3-dimensional ndarray indexed by spherical index, radial index, and time index

Returns

numpy.ndarray: 3-dimensional numpy array containing the spherical coefficients

Notes

The spherical index serializes all pairs of (l, m). The index for (l, m) is calculated as: l*(l+1)/2 + m.

setMatrix(self: pyEXP.coefs.SphCoefs, time: float, mat: numpy.ndarray[numpy.complex128[m, n]]) → None

Enter and/or rewrite the coefficient matrix at the provided time

Parameters

timefloat: snapshot time corresponding to the the coefficient matrix mat : numpy.ndarray

the new coefficient array.

Returns

None

class pyEXP.coefs.SphFldCoefs

Container for spherical field coefficients

getAllCoefs(self: pyEXP.coefs.SphFldCoefs) → numpy.ndarray[numpy.complex128]

Provide a 4-dimensional ndarray indexed by channel index, spherical index, radial index, and time index

Returns

numpy.ndarray: 4-dimensional numpy array containing the spherical coefficients

Notes

The spherical index serializes all pairs of (l, m) where l, m are the aximuthal indices. The index for (l, m) pair is calculated as: l*(l+1)/2 + m

setMatrix(self: pyEXP.coefs.SphFldCoefs, time: float, mat: numpy.ndarray[numpy.complex128]) → None

Enter and/or rewrite the coefficient tensor at the provided time

Parameters

timefloat: snapshot time corresponding to the the coefficient matrix
matnumpy.ndarray: the new coefficient array.

Returns

None

class pyEXP.coefs.TableData

Container for simple data tables with multiple columns

getAllCoefs(self: pyEXP.coefs.TableData) → numpy.ndarray[numpy.float64[m, n]]

Return a 2-dimensional ndarray indexed by column and time

Returns

numpy.ndarray: 2-dimensional numpy array containing the data table

class pyEXP.coefs.TrajectoryData

Container for trajectory/orbit data

getAllCoefs(self: pyEXP.coefs.TrajectoryData) → numpy.ndarray[numpy.float64]

Return a 3-dimensional ndarray indexed by column and time

Returns

numpy.ndarray: 2-dimensional numpy array containing the data table

FieldGenerator class bindings

FieldGenerator

This class computes surfaces and volumes for visualizing the physical quantities implied by your basis and associated coefficients.

The generator is constructed by passing a vector of desired times that must be in the coefficient object and list of lower bounds, a list of upper bounds, and a list of knots per dimension. These lists all have rank 3 for (x, y, z). For a two-dimensional surface, one of the knot array values must be zero. The member functions lines, slices, an arbitrary point-set, and volumes, called with the basis and coefficient objects, return a numpy.ndarray containing the field evaluations. Each of these functions returns a dictionary of times to a dictionary of field names to numpy.ndarrays at each time. There are also members which will write these generated fields to files. The linear probe members, ‘lines’ and ‘file_lines’, evaluate ‘num’ field points along a user-specified segment between the 3d points ‘beg’ and ‘end’. See help(pyEXP.basis) and help(pyEXP.coefs) for info on the basis and coefficient objects.

Data packing

All slices and volumes are returned as numpy.ndarrays in row-major order. That is, the first index is x, the second is y, and the third is z. The ranks are specified by the ‘gridsize’ array with (nx, ny, nz) as input to the FieldGenerator constructor. A point mesh is rows of (x, y, z) Cartesian coordinates. The output field values returned by points is in the same order as the input array.

Coordinate systems

The FieldGenerator class supports spherical, cylindrical, and Cartesian for force field components. These are selected by your basis instance, consistent with the natural coordinate system for that basis. You may change the default coorindate system for a basis using the ‘setFieldType()’ member function.

Field names

The data fields are as follows: * ‘dens’ the total density * ‘dens m>0’ the non-axisymmetric component of the density * ‘dens m=0’ the axisymmetric component of the density * ‘rad force’ the radial force * ‘mer force’ the meridional force * ‘azi force’ the azimuthal force * ‘potl’ the total potential * ‘potl m>0’ the non-axisymmetric component of the potential * ‘potl m=0’ the axisymmetric component of the potential

For spherical coordinates (coord=”Spherical”): * ‘rad force’ the radial force * ‘mer force’ the meridional force * ‘azi force’ the azimuthal force

For cylindrical coordinates (coord=”Cylindrical”): * ‘rad force’ the radial force * ‘ver force’ the meridional force * ‘azi force’ the azimuthal force

For Cartesian coordinates (coord=”Cartesian”): * ‘x force’ the radial force * ‘y force’ the meridional force * ‘z force’ the azimuthal force

Notes

Note that the ‘dens’ field is the sum of the ‘dens m=0’ and ‘dens m>0’ fields and, similarly, the ‘potl’ field is the sum of ‘potl m=0’ and ‘potl m>0’ fields. These redundant entries are provided for convenience and conceptual clarity. For spherical bases, the ‘m’ index should be interpreted as the ‘l’ index.

Multivariate Singular Spectrum Analysis (MSSA) class bindings

expMSSA

The expMSSA class analyzes and separates the temporal patterns in your coefficients and auxiliary data. Instances can return coefficient sets for one or more individual patterns that may be used with the FieldGenerator for visualization

Coefficient update and memory usage

A key feature of MSSA inference is the field generation from the coefficients reconstructed from dominant eigenvalues. After identifying interesting PCs and associated eigenvalues, this is done in ‘expMSSA’ in two steps:

A call to ‘reconstruct(list)’ with a list of eigenvalue indices reconstruct the data series and saves those series to working (possibly detrended) vectors

2 ‘getReconstructed()’ returns a dictionary of name strings that
point to coefficient objects (Coefs). These may be used for diagnostics (e.g. visualizing fields).

The split between the the ‘reconstruct()’ and the ‘getReconstructed() members allows one to examine the detrended series with the selected eigenvalues before updating the coefficient data. The ‘reconstruct() call creates the internal working channel data, leaving the initial coefficient set untouched. The ‘getReconstructed()’ call replaces the coefficient data in the Coefs object passed to expMSSA, as well as returning a dictionary of name strings and coefficient containers. Subsequent calls to ‘getReconstructed’ will overwrite the previous updates. In practice, the coefficients object that you made in Python will always contain the latest updates from the last ‘getReconstructed’ call. This allows you to update channels incrementally and the updated copy is available to Python as the originally passed reference. This strategy prevents your stack from growing very large by efficiently reusing previously allocated storage. If you want to save a copy of the original coefficients, use the ‘deepcopy()’ member for Coefs before the reconstruction. The original ‘factory’ call and each call to ‘deepcopy()’ will make a new set that will be freed when the reference to the coefficient instance disappears. This allows some explicit control over your memory use.

Configuration parameters

expMSSA has a number of internal configuration parameters, listed below. Most people will not need them, but they do exist should you really want to dig in.

To change the default expMSSA, pass a YAML database of parameters to the expMSSA constructor. This is a convenient way to have control over some fine-tuning parameters or experimental features without a a change to the API. They are all optional but a few of these will be useful. They are all YAML key-value pairs. The first group below are simple boolean toggles; their presence turns on a feature independent of the value. The Python script samples show how to construct these simple YAML configurations on the fly. A simple example is also given below. The boolean parameters are listed below by my guess of their usefulness to most people:

verbose: false Report on internal progress for debugging only noMean: false If true, do not subtract the mean when

reading in channels. Used only for totPow detrending method.

writeCov: true Write the covariance matrix to a file for
diagnostics since this is not directly available from the interface. Used only if Traj: false.

Jacobi: true Use the Jacobi SVD rather than the Random
approximation algorithm from Halko, Martinsson, and Tropp (RedSVD). This is quite accurate but _very_ slow

BDCSVD: true Use the Bidiagonal Divide and Conquer SVD
rather than RedSVD; this is faster and more accurate than the default RedSVD but slower.

Traj: true Perform the SVD of the trajectory matrix
rather than the more computationally intensive but more stable SVD of the covariance matrix. Set to false to get standard covariance SVD. In practice, ‘Traj: true’ is sufficiently accurate for all PC orders that are typically found to have dynamical signal.

rank: 100 The default rank for the randomized matrix SVD.
The default value will give decent accuracy with small computational overhead and will be a good choice for most applications.

RedSym: true Use the randomized symmetric eigenvalue solver
RedSym rather rather than RedSVD for the co- variance matrix SVD (Traj: false). The main use for this is checking the accuracy of the default randomized matrix methods.

allchan: true Perform k-means clustering analysis using all
channels simultaneously

distance: true Compute w-correlation matrix PNG images using
w-distance rather than correlation

flip: true Exchanges the x-y axes in PNG plots

power: false Compute and output Fourier power totVar: false Detrend according to the total variance

in all channel

totPow: false Detrend according to the total power in
all channels

The following parameters take values, defaults are given in ()

evtol: double(0.01) Truncate by the given cumulative p-value in
chatty mode

output: str(exp_mssa) Prefix name for output files

The ‘output’ value is only used if ‘writeFiles’ is specified, too. A simple YAML configuration for expMSSA might look like this: — BDCSVD: true evtol: 0.01 …

If you find that you are using these YAML parameters often, let me know and I can promote them to the main API

Grouping

MSSA often spreads the same signal between more than one PC. We try to group eigenvalue-PC pairs that seem to have similar temporatl behavior in the final stage before reconstruction. The w-correlation matrix, the k-means analyses and similarities of the PC frequency spectra are are ways of diagnosing groups of related components. You may group using ‘reconstruct(list)’, where ‘list’ is a Python list of component indices in eigenvalue-PC index order. An empty argument implies no that the will be no PCs in the reconstruction, only the mean values. The two expMSSA steps will often need to be used iteratively. First, to learn about the primary signals in your data and catagorize the PCs into groups. Second, a new analysis that provides separation of your main distinct signals into individual reconstructions.

Save/Restore

MSSA analyses can be computationally intensive so we include a way to serialize the state to an HDF5 file. The saveState(prefix) member takes a prefix file parameter and creates the file <prefix>_mssa.h5. You may restore the state by recreating your expMSSA instance with the identical input data and keys and then using restoreState(prefix) to read the MSSA analysis from the HDF5 file. The restore step will check that your data has the same dimension, same parameters, and key list so it should be pretty hard to fool, but not impossible.

Computational notes

Some of the linear algebra operations and the MSSA reconstruction can parallelize itself using OpenMP threads. If your compilter has OpenMP, please enable it (i.e. -fopenmp when compiling in GNU or and set the appropriate environment (e.g. ‘export OMP_NUM_THREADS=X’ for Bash). This will be especially useful for the ‘reconstruct()’ step. The initial MSSA analysis can also be memory intensive, depending on the size of the time series and the number of channels. One possible strategy is to run the analyses on a larger memory compute node and save the results using the ‘saveState()’ member to an HDF5 file and reread those files on a local machine using the ‘restoreState()’ member.

Extended Dynamical Mode Decomposition (EDMD) class bindings

Koopman

This Koopman class implements EDMD, using the exact DMD algorith from Tu et al. (2014) using the EXP coefficients and your optional auxiliary data. Just as in expMSSA, instances return coefficient sets for one or more individual EMD modes. These may, in turn, be used with the FieldGenerator for visualization

Coefficient update and memory usage

A key feature of Koopman operator theory is the identification of the underlying dynamical structure in the measurement of state space. In this case, the measurement of phase space is the collection of scalars that represent the density-potential fields. There is the possibility and maybe even the likelihood that the mismatch of the fields to the true nature of the dynamics will generate spurious models. This is true with mSSA as well. In fact, mSSA is also an approximation to the Koopman operator, so the philosophy is similar. My experience to date suggests that mSSA has much better performance and separation of signals. If you try this and have good success, please do contact me. As as in in mSSA, this class allows you to gain insight from a set of coefficients reconstructed from dominant EDMD modes. After identifying interesting modes and associated eigenvalues, this is done in ‘Koopman’ in two steps:

A call to ‘reconstruct(list)’ with a list of eigenvalue indices reconstruct the data series and saves those series to working (possibly detrended) vectors

2 ‘getReconstructed()’ returns a dictionary of name strings that
point to coefficient objects (Coefs). These may be used for diagnostics (e.g. visualizing fields).

This structure exactly parallels the implementation in ‘expMSSA’ The split between the the ‘reconstruct()’ and the ‘getReconstructed() members allows one to examine the detrended series with the selected eigenvalues before updating the coefficient data. The ‘reconstruct() call creates the internal working channel data, leaving the initial coefficient set untouched. The ‘getReconstructed()’ call replaces the coefficient data in the Coefs object passed to Koopman, as well as returning a dictionary of name strings and coefficient containers. Subsequent calls to ‘getReconstructed’ will overwrite the previous updates. In practice, the coefficients object that you made in Python will always contain the latest updates from the last ‘getReconstructed’ call. This allows you to update channels incrementally and the updated copy is available to Python as the originally passed reference. This strategy prevents your stack from growing very large by efficiently reusing previously allocated storage. If you want to save a copy of the original coefficients, use the ‘deepcopy()’ member for Coefs before the reconstruction. The original ‘factory’ call and each call to ‘deepcopy()’ will make a new set that will be freed when the reference to the coefficient instance disappears. This allows some explicit control over your memory use.

Configuration parameters

Koopman has a number of internal configuration parameters, listed below. Most people will not need them, but they do exist should you really want to dig in. These are a subset of the same parameter available in expMSSA.

To change the default values, pass a YAML database of parameters to the Koopman constructor. This is a convenient way to have control over some fine-tuning parameters or experimental features without a a change to the API. They are all optional but a few of these will be useful. They are all YAML key-value pairs. The first group below are simple boolean toggles; their presence turns on a feature independent of the value. The Python script samples show how to construct these simple YAML configurations on the fly. A simple example is also given below. The boolean parameters are listed below by my guess of their usefulness to most people:

verbose: false Whether there is report or not Jacobi: true Use the Jacobi SVD rather than the Random

approximation algorithm from Halko, Martinsson, and Tropp (RedSVD). This is quite accurate but _very_ slow

BDCSVD: true Use the Binary Divide and Conquer SVD rather
rather than RedSVD; this is faster and more accurate than the default RedSVD but slower

project: true Use the classic DMD projected modes rather than
the exact DMD modes, as described by Tu et al. 2014.

power: true Write partial power contributions into a file in
a ascii table format if set to ‘true’. Default is ‘false’

The following parameters take values, defaults are given in ()

output: sting Prefix name for output files. The default is
‘exp_edmd’.

The ‘output’ value is used by ‘getContributions()’ and ‘channelDFT()’ if the ‘power’ options is set. A simple YAML configuration for Koopman might look like this: — BDCSVD: true project: true …

If you find that you are using these YAML parameters often, let me know and I can promote them to the main API

Save/Restore

Koopman analyses can be computationally intensive so we include a way to serialize the state to an HDF5 file. The saveState(prefix) member takes a prefix file parameter and creates the file <prefix>_edmd.h5. You may restore the state by recreating your Koopman instance with the identical input data and keys and then using restoreState(prefix) to read the Koopman analysis from the HDF5 file. The restore step will check that your data has the same dimension, same parameters, and key list so it should be pretty hard to fool, but not impossible. The computation is much less expensive, generally, than MSSA so saving the analysis is not strictly necessary, even on a laptop.

Computational notes

Some of the linear algebra operations can parallelize itself using OpenMP threads. If your compilter has OpenMP, please enable it (i.e. -fopenmp when compiling in GNU or and set the appropriate environment (e.g. ‘export OMP_NUM_THREADS=X’ for Bash). The initial Koopman analysis can also be memory intensive, depending on the size of the time series and the number of channels. Although it requires much less memory than the equivalent mSSA analysis. If memory is a problem, try running the analyses on a larger memory compute node and save using the ‘saveState()’ member to an HDF5 file and reread those files on a local machine using the ‘restoreState()’ member.

Utility class bindings

This module provides routines for BFE tasks that do not naturally fit into the main categories. The current list of utilities is:

Report the current EXP version, GIT branch and commit

Compute the center of the particle distribution from its center of mass. Very fast but easily biased.

Compute the mean density weighted center of the particle distribu- tion from KD density estimator at each particle position. Very is very slow. One can change the stride to decrease the sample size to speed this up.

Note on centers: the EXP n-body code does this automatically. The density weighted center is an alternative for snapshots without center estimates. Only use COM if know your simulation remains close to bisymmetric.

Apply a user-defined Python function to all particles in a phase- space reader. This may be used to do calculations using all or a user-determined subset of snapshot particles. The functor has no return type; it is up to the user to put accumulated values in the scope. The functor needs the arguments of mass, position vector, velocity vector, and index. For example, the following Python code computes the center of mass: #————————————————————— # Variables in the scope of myFunctor # totalMass = 0.0 centerOfMass = [0.0, 0.0, 0.0] # # Definition of the functor # def myFunctor(m, pos, vel, index):

global totalMass, centerOfMass totalMass += m for i in range(3): centerOfMass[i] += m*pos[i]

# # Now iterate through the particles provided by a reader instance # pyEXP.util.particleIterator(reader, myFunctor) # # Print the COM # for i in range(3): centerOfMass[i] /= totalMass # #—————————————————————

pyEXP.util.Version() → tuple[int, int, int]

Return the version.

This is the same version information reported by the EXP N-body code.

pyEXP.util.getCenterOfMass(reader: pyEXP.read.ParticleReader) → list[float]

Compute the center of mass for the particle component

Parameters

readerParticleReader: the particle-reader class instance

Returns

list(float): Computed center

pyEXP.util.getDensityCenter(reader: pyEXP.read.ParticleReader, stride: int = 1, Nsort: int = 0, Ndens: int = 32) → list[float]

Compute the center of the particle component

This implementation uses the density weighted position using KD N-nearest neighbor estimator.

Parameters

readerParticleReader: the particle-reader class instance
strideint, default=1: stride >1 will generate a subsample of every nth particle over a random permutation
Ndensint, default=32: number of particles per sample ball (32 is a good compromise between accuracy and runtime; 16 is okay if you are trying to shave off runtime.
Nsortint, default=0: Nsort >0 keeps the particles of the Nsort densest samples

Returns

list(float): Computed center

pyEXP.util.getVersionInfo() → None

Report on the version and git commit.

This is the same version information reported by the EXP N-body code.

pyEXP.util.particleIterator(reader: pyEXP.read.ParticleReader, functor: Callable[[float, list[float], list[float], int], None]) → None

Apply a user-defined functor to every particle in phase space

Parameters

readerParticleReader: the particle-reader class instance
functor :: the callback function to compute a phase-space quantity

Returns

None

Notes

The callback function must have the signature:

void(float, list, list, int)

where the first argument is mass, the second is position, the third is velocity, and the fourth is index. Not all values need to be used in the function, of course.