An example of creating coefficients using pyEXP

This example uses a prerun simulation whose output stored in the Tutorials/Data directory. This can be recreated using EXP and the example configuration and input files from the EXP-examples/Nbody simulation. But this notebook could be adapted for any simulation you like.

We begin by importing pyEXP and friends and setting the working directory.

import os
import yaml
import pyEXP

# As described in the README, we assume that you have started Jupyter in the `Tutorials` directory and have provided all of the
# necessary data in the `Data` subdirectory.  We move to that directory as a first step:
#
os.chdir('../Data')

Create the basis

We’ll only do the halo coefficients in this simple example. The cylindrical coefficients would procede similarly. See the viewing a basis notebook for an example of creating the cylindrical basis.

# Get the basis config.  This loads a YAML stanza specifies all of the input parameters
# necessary to build the basis.
#
yaml_config = ""
with open('basis.yaml') as f:
    config = yaml.load(f, Loader=yaml.FullLoader)
    yaml_config = yaml.dump(config)

# Alternatively, you could construct this on the fly in Python, e.g.
bconfig = """
---
id: sphereSL
parameters :
  numr: 2000
  rmin: 0.0001
  rmax: 1.95
  Lmax: 4
  nmax: 10
  rmapping: 0.0667
  modelname: SLGridSph.model
  cachename: sphereSL.cache
...
"""
print('-'*60)
print('Read from file')
print('-'*60)
print(yaml_config)
print('-'*60)
print('Constructed from string')
print('-'*60)
print(bconfig)
print('-'*60)

# Construct the basis instance
#
basis   = pyEXP.basis.Basis.factory(yaml_config)

------------------------------------------------------------
Read from file
------------------------------------------------------------
id: sphereSL
parameters:
  Lmax: 4
  cachename: sphereSL.model
  modelname: SLGridSph.model
  nmax: 10
  numr: 2000
  rmapping: 0.0667
  rmax: 1.95
  rmin: 0.0001

------------------------------------------------------------
Constructed from string
------------------------------------------------------------

---
id: sphereSL
parameters :
  numr: 2000
  rmin: 0.0001
  rmax: 1.95
  Lmax: 4
  nmax: 10
  rmapping: 0.0667
  modelname: SLGridSph.model
  cachename: sphereSL.cache
...

------------------------------------------------------------
---- SLGridSph::ReadH5Cache: successfully read basis cache <sphereSL.model>
---- Spherical::orthoTest: worst=0.00016446

When an adaptive basis such as sphereSL or Cylindrical is constrructed, the orthogonality of the biorthogonal basis is automatically computed as a sanity check on the model file and input parameters. The condition reported by orthoTest is the absolute value of the inner product minus one or the inner product itself, for diagonal and off-diagonal matrix elements, respectively. In this case, you can see that the worst inner product is ~0.0001 or on part in \(10^4\).

Creating a particle reader

Now that we have a basis, we can use it to create coefficients from the particle snapshots. pyEXP uses a ParticleReader object for that.

The first step is to hand off the files that comprise a snapshot for every time slice. The ParticleReader provides a helper function for that. There are two helper functions: parseFileList and parseStringList. The first reads a list from a file and the second takes a list. Otherwise they are the same. The file names in the list are assumed to end with a snapshot index and an optional part index. For example, if you have single files per snapshot, the list might look like: myrun.00000, myrun.00001, etc. If you have multiple files per snapshot, they will look something like myrun.00000_0001, myrun.00000_0002, myrun.00001_0000, myrun.00001_0001, etc.

Here is the call for a file:

# Construct batches of files the particle reader.  One could use the
# parseStringList to create batches from a vector/list of files.  NB:
# a std::vector in C++ becomes a Python.list and vice versa
#
batches = pyEXP.read.ParticleReader.parseFileList('file.list', '')

We now iterate the batches created by the file parser to create the coefficients. For each batch we create a new reader and pass the reader to the basis instance. The basis.createFromReader member creates and returns the coefficients. The coefficients are added to a coefficient container called coefs. Note: on the first call coefs=None so a new container is created on the first time through.

# This will contain the coefficient container, need to start with a
# null instance to trigger construction
#
coefs   = None

for group in batches:

    print("file group is", group)

    # Make the reader for the desired type.  One could probably try to
    # do this by inspection but that's another project.
    #
    reader = pyEXP.read.ParticleReader.createReader('PSPout', group, 0, False);

    # Print the type list
    #
    print('The component names are:', reader.GetTypes())

    compname = 'dark halo'
    reader.SelectType(compname)
    print('Selected', compname)

    print('Call createFromReader at Time', reader.CurrentTime(), 'for', reader.CurrentNumber(), 'particles')

    coef = basis.createFromReader(reader)
    print("Created coef")

    # We need this idiom here because None is not mapping to a
    # null pointer in pybind11.  There is probably a way to do this.
    # Suggestions anyone?
    #
    #                          This is optional---+
    #                                             |
    if coefs is None:           #                 v
        coefs = pyEXP.coefs.Coefs.makecoefs(coef, compname)
    else:
        coefs.add(coef)

    print('Added coef')
    print('-'*60)

print('\nCompleted the file group list\n')

print('The coefficient time list is', coefs.Times())

file group is ['OUT.run0.00000']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.0 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00001']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.009999999999999764 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00002']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.019999999999998665 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00003']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.029999999999997563 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00004']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.04000000000000035 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00005']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.05000000000000369 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00006']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.06000000000000703 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00007']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.0700000000000037 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00008']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.07999999999999816 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00009']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.08999999999999261 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00010']
The component names are: ['dark halo', 'star disk']
Selected dark halo
 for 100000 particles at Time 0.09999999999998707
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00011']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.10999999999998153 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00012']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.11999999999997599 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00013']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.12999999999997933 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00014']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.13999999999999155 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00015']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.15000000000000377 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00016']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.160000000000016 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00017']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.1700000000000282 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00018']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.18000000000004043 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00019']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.19000000000005265 for 100000 particles
Created coef
Added coef
------------------------------------------------------------
file group is ['OUT.run0.00020']
The component names are: ['dark halo', 'star disk']
Selected dark halo
Call createFromReader at Time 0.20000000000006488 for 100000 particles
Created coef
Added coef
------------------------------------------------------------

Completed the file group list

The coefficient time list is [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2]

The reader gives a verbose status update for each snapshot as it is read. You can see the component names, simulation time, and particles numbers per snapshot as the files are read. For large simulations, coefficient construction can take awhile. Check out the following two scripts in the Utitilies directory:

1. make coefficients MPI.py

2. make coefficients native MPI.py

The first one is an example for Gadget snapshots and the second is an example for EXP snapshots.

Using a FieldGenerator

Now that we have our new coefficients, we can use the FieldGenerator to view the BFE representation of the underlying fields. Here is an example:

# Now try some slices for rendering
#

times = coefs.Times()[0:3]
pmin  = [-1.0, -1.0, 0.0]
pmax  = [ 1.0,  1.0, 0.0]
grid  = [  40,   40,   0]

fields = pyEXP.field.FieldGenerator(times, pmin, pmax, grid)

surfaces = fields.slices(basis, coefs)

print("We now have the following [time field] pairs")
for v in surfaces:
    print('-'*40)
    for u in surfaces[v]:
        print("{:8.4f}  {}".format(v, u))

print("\nHere is the first one:")
for v in surfaces:
    for u in surfaces[v]:
        print('-'*40)
        print('----', u)
        print('-'*40)
        print(surfaces[v][u])
    break

We now have the following [time field] pairs
----------------------------------------
  0.0100  azi force
  0.0100  dens
  0.0100  dens m=0
  0.0100  dens m>0
  0.0100  mer force
  0.0100  potl
  0.0100  potl m=0
  0.0100  potl m>0
  0.0100  rad force
----------------------------------------
  0.0200  azi force
  0.0200  dens
  0.0200  dens m=0
  0.0200  dens m>0
  0.0200  mer force
  0.0200  potl
  0.0200  potl m=0
  0.0200  potl m>0
  0.0200  rad force
----------------------------------------
  0.0300  azi force
  0.0300  dens
  0.0300  dens m=0
  0.0300  dens m>0
  0.0300  mer force
  0.0300  potl
  0.0300  potl m=0
  0.0300  potl m>0
  0.0300  rad force

Here is the first one:
----------------------------------------
---- azi force
----------------------------------------
[[-0.00087384 -0.00098027 -0.00113208 ...  0.00053854  0.00058888
   0.00065692]
 [-0.00104563 -0.00117976 -0.00136719 ...  0.00044353  0.00050006
   0.00057894]
 [-0.00120132 -0.0013625  -0.00158698 ...  0.00027801  0.00034734
   0.00044279]
 ...
 [-0.00086897 -0.00120733 -0.00153372 ...  0.00364213  0.00354581
   0.00338669]
 [-0.00081985 -0.00113707 -0.00144417 ...  0.00372749  0.00363169
   0.00347801]
 [-0.00079396 -0.00109152 -0.00138081 ...  0.00380656  0.00370932
   0.00355781]]
----------------------------------------
---- dens
----------------------------------------
[[0.01483475 0.01596043 0.01714892 ... 0.01672023 0.01557306 0.01447723]
 [0.01600303 0.01726678 0.01860979 ... 0.01800252 0.01671672 0.01549673]
 [0.01722907 0.01864478 0.02016317 ... 0.01938812 0.01793789 0.01657638]
 ...
 [0.01678794 0.01815355 0.01959459 ... 0.02029726 0.01882576 0.01744723]
 [0.01564873 0.01686448 0.01813334 ... 0.0186915  0.01740002 0.01617756]
 [0.01458028 0.01566568 0.0167905  ... 0.01718854 0.01605308 0.01497049]]
----------------------------------------
---- dens m=0
----------------------------------------
[[0.01438968 0.0155077  0.01669507 ... 0.01669507 0.0155077  0.01438968]
 [0.0155077  0.01676509 0.01810998 ... 0.01810998 0.01676509 0.0155077 ]
 [0.01669507 0.01810998 0.01963967 ... 0.01963967 0.01810998 0.01669507]
 ...
 [0.01669507 0.01810998 0.01963967 ... 0.01963967 0.01810998 0.01669507]
 [0.0155077  0.01676509 0.01810998 ... 0.01810998 0.01676509 0.0155077 ]
 [0.01438968 0.0155077  0.01669507 ... 0.01669507 0.0155077  0.01438968]]
----------------------------------------
---- dens m>0
----------------------------------------
[[ 4.45073500e-04  4.52731590e-04  4.53846878e-04 ...  2.51552956e-05
   6.53618481e-05  8.75497790e-05]
 [ 4.95328044e-04  5.01695555e-04  4.99815797e-04 ... -1.07452950e-04
  -4.83680560e-05 -1.09646717e-05]
 [ 5.34000923e-04  5.34801453e-04  5.23498224e-04 ... -2.51551945e-04
  -1.72089174e-04 -1.18696655e-04]
 ...
 [ 9.28678128e-05  4.35745387e-05 -4.50829830e-05 ...  6.57595694e-04
   7.15783623e-04  7.52157473e-04]
 [ 1.41032098e-04  9.93959984e-05  2.33634619e-05 ...  5.81527303e-04
   6.34925615e-04  6.69857895e-04]
 [ 1.90603168e-04  1.57983566e-04  9.54256029e-05 ...  4.93466156e-04
   5.45385177e-04  5.80809254e-04]]
----------------------------------------
---- mer force
----------------------------------------
[[ 0.00116107  0.0012984   0.00143422 ...  0.00141841  0.00133688
   0.00125844]
 [ 0.0010585   0.00120293  0.00134758 ...  0.00156861  0.00148209
   0.00139784]
 [ 0.00094972  0.00110494  0.00126318 ...  0.00172921  0.00163881
   0.00154933]
 ...
 [-0.00298936 -0.003173   -0.00336324 ...  0.00119402  0.00105672
   0.00093735]
 [-0.00276366 -0.00292988 -0.0031003  ...  0.00081658  0.00071328
   0.00062615]
 [-0.00254213 -0.00269264 -0.00284517 ...  0.0004661   0.0003963
   0.00034138]]
----------------------------------------
---- potl
----------------------------------------
[[-0.9169819  -0.9383662  -0.95999885 ... -0.960537   -0.9389293
  -0.917559  ]
 [-0.9383129  -0.9611787  -0.9843996  ... -0.984946   -0.96175075
  -0.93889844]
 [-0.95987254 -0.98432374 -1.0092399  ... -1.009817   -0.98492396
  -0.960483  ]
 ...
 [-0.9592484  -0.98363    -1.0084537  ... -1.0091766  -0.984431
  -0.9601257 ]
 [-0.93764687 -0.9604352  -0.98355633 ... -0.98422873 -0.9611766
  -0.9384567 ]
 [-0.91629386 -0.9375966  -0.9591268  ... -0.9597348  -0.9382673
  -0.9170273 ]]
----------------------------------------
---- potl m=0
----------------------------------------
[[-0.9166649 -0.9380169 -0.9596112 ... -0.9596112 -0.9380169 -0.9166649]
 [-0.9380169 -0.9608532 -0.9840379 ... -0.9840379 -0.9608532 -0.9380169]
 [-0.9596112 -0.9840379 -1.0089229 ... -1.0089229 -0.9840379 -0.9596112]
 ...
 [-0.9596112 -0.9840379 -1.0089229 ... -1.0089229 -0.9840379 -0.9596112]
 [-0.9380169 -0.9608532 -0.9840379 ... -0.9840379 -0.9608532 -0.9380169]
 [-0.9166649 -0.9380169 -0.9596112 ... -0.9596112 -0.9380169 -0.9166649]]
----------------------------------------
---- potl m>0
----------------------------------------
[[-0.00031698 -0.00034926 -0.00038765 ... -0.00092586 -0.00091239
  -0.00089413]
 [-0.00029595 -0.00032549 -0.00036175 ... -0.00090816 -0.00089753
  -0.00088152]
 [-0.00026139 -0.00028588 -0.00031703 ... -0.00089411 -0.00088608
  -0.00087182]
 ...
 [ 0.00036277  0.00040785  0.00046916 ... -0.00025367 -0.00039316
  -0.00051454]
 [ 0.00037004  0.00041801  0.00048153 ... -0.00019086 -0.00032336
  -0.00043978]
 [ 0.00037099  0.00042032  0.00048442 ... -0.00012365 -0.00025039
  -0.0003624 ]]
----------------------------------------
---- rad force
----------------------------------------
[[-0.58504117 -0.60900503 -0.63369316 ... -0.63327456 -0.6087623
  -0.5849374 ]
 [-0.6087843  -0.6348392  -0.66203845 ... -0.66186184 -0.6348188
  -0.60887575]
 [-0.6331737  -0.661728   -0.691743   ... -0.6919265  -0.66203153
  -0.6335524 ]
 ...
 [-0.6329953  -0.6614362  -0.6912934  ... -0.69079643 -0.66104454
  -0.632722  ]
 [-0.60824454 -0.6341394  -0.66113216 ... -0.66103613 -0.6340929
  -0.6082666 ]
 [-0.58422893 -0.6079955  -0.6324451  ... -0.632654   -0.6082124
  -0.5844728 ]]

These could be make into images and so forth. We’ll do this in another example notebook.

Saving the coefficients

At this point, it makes sense to save the coefficients that you have just created. This is sone with the following call:

# Remove any old file with the same name
try:
    os.remove('test_coefs')
except:
    print('Nothing to delete')
coefs.WriteH5Coefs('test_coefs')
print('Coefficient file saved')

Nothing to delete
Coefficient file saved

We now have a EXP HDF5 coefficient file called test_coefs. You can view the contents of test_coefs directly using the h5dump command supplied in your HDF5 installion:

h5dump test_coefs | less

As an example of manipulating the newly made coefficients in pyEXP, let’s try reading the newly created file into another coefficient container, coefs2. The container has a member function called CompareStanzas which will check on the contents. Let’s do it.

# Now try reading it in
#
coefs2 = pyEXP.coefs.Coefs.factory('test_coefs')
print("Type is", coefs2.getGeometry())

# Now compare with the original
#
coefs2.CompareStanzas(coefs)

Type is sphere
Times are the same, now checking parameters at each time
Parameters are the same, now checking coefficients

True

This member function will print differences. No differenced should be printed, of course.

---

Where do you want to go next?

• Return to Table of Contents

• Go to next notebook