Add some units documentation

index.rst

@@ -89,6 +89,7 @@ EXP concepts
    topics/bfetheory
    topics/codeintro
    topics/timesseries
+   topics/units
    topics/multistep
    topics/centering
    topics/yamlconfig
@@ -106,6 +107,9 @@ mathematics used in EXP.
 :doc:`topics/timesseries`
     Descriptions for methods to analyze BFE time series
 
+:doc:`topics/units`
+    Physical units in EXP 
+
 :doc:`topics/multistep`
     A quick review of EXP's ODE solver and multi time step ladder
 

intro/overview.rst

@@ -96,6 +96,11 @@ Here's the code for computing coefficients::
      #
      coef = basis.createFromArray(data[:,0], data[:,1:4], time=3.0)
 
+     # Add your simulation units (we assume a unit-free simulation here)
+     #
+     coef.setUnits([('mass', 'none', 1.0), ('length', 'none', 1.0),
+     ('time', 'none', 1.0), ('G', 'none', 1.0)])
+
      # Make an HDF5 file
      #
      coefs = pyEXP.coefs.Coefs.makecoefs(coef)
@@ -144,6 +149,11 @@ More detailed information on YAML and config parameters is available
 in the :ref:`What is YAML?<yamlconfig>` and :ref:`How to visualize the
 BFE bases used to make your coefficients<visualizing-bases>` pages.
 
+The unit system in EXP is described in more detail in :ref:`Units in
+EXP and pyEXP<units>`.  For this example, we assume a unit-free
+simulation but it's easy to add any unit system that is convenient for
+you.
+
 pyEXP is then ready to make the coefficients from your phase-space
 data.  This example assumes that the mass and positions of your
 particles are in columns 1, 2, 3, 4 of the file and that the positions

topics/units.rst

@@ -0,0 +1,291 @@
+.. role:: python(code)
+       :language: python
+       :class: highlight
+
+.. _units:
+
+Units in EXP and pyEXP
+======================
+
+Overview
+--------
+
+The EXP N-body code assumes the gravitational constant is unity and
+that mass, position and velocity units can have any value defined by
+the user.  In other words, the EXP N-body code does not enforce any
+particular unit set consistent with ``G=1``.
+
+pyEXP will also accept mass, position, time, and velocity with any
+unit set defined by the imported simulation as well as an arbitrary
+value of the gravitational constant. Explicit units types and the
+gravitational constant ``G`` may be set at construction time as we
+will describe below in :ref:`unit-schema`.  pyEXP **requires** the
+user to set units explicitly when coefficient sets are
+constructed. There is no default. See :ref:`intro-pyEXP-tutorial` for
+an end-to-end example of coefficient computation.
+
+All unit information is written into the EXP coefficient files as HDF5
+attributes.  Also see :ref:`hdf5-unit-attributes` for details.  You
+may add units to previously computed coefficient files using the
+script described in :ref:`units_old`.
+
+.. _unit-schema:
+
+The EXP unit schema
+-------------------
+
+EXP requires one of the following two sequences of 4 unit types:
+
+1. Mass, Length, Time, gravitational constant (G)
+
+2. Mass, Length, Velocity, gravitational constant (G)
+
+
+Other combinations are possible in principle, such as Mass, Length,
+Velocity and Time.  However, this would require an automatic deduction
+of the gravitational constant. EXP does not current know about
+physical unit conversion.  This feature may be added in the future.
+
+Each separate unit is defined as a ``tuple`` which takes the form::
+
+  ('unit type', 'unit name', <float value>)
+
+The type and name strings are checked against the allowed values as follows:
+
+- The ``unit type`` is one of 'length', 'mass', 'time', 'velocity' or
+  'G'. The internal parser will accept a variety of aliases for these
+  such as 'Length', 'Len', 'len', 'l', 'L' for 'length'; 'Mass', 'm',
+  'M' for 'mass'; 'Time', 't', 'T' for 'time', 'vel',
+  'Vel', 'Velocity', 'v', 'V' for 'velocity'; and 'Grav', 'grav',
+  'grav_constant', 'Grav_constant', 'gravitational_constant',
+  'Gravitational_constant' for 'G'.
+
+- The ``unit name`` is one of the usual unit names for each of the
+  ``unit type``.  The allowed list is a subset of the standard
+  `astropy` strings.  For example, 'pc' or 'kpc' for 'length'.  There
+  is also a 'none' type which implies no assigned physical units.
+
+- The floating value defines the number of each unit for the type.
+  The value can be any valid floating-point number.
+
+The allowed types and names may be queried interactively in pyEXP
+using the :python:`getAllowedUnitTypes()` and
+:python:`getAllowedUnitNames(type)`, see :ref:`units-interface`.  For
+development reference, these allowed strings are defined in
+``expui/UnitValidator.cc`` in the EXP source tree.
+
+As an example, a Gadget user might define mass units as::
+
+  ('mass', 'Msun', 1.0e10)
+
+The full units specification is a list of tuples that includes one of
+the four sequences.  In C++, the list is a :cpp:any:`std::vector<>` of
+tuples.
+
+As an example, all native EXP simulations have the following unit
+lists:
+
+   .. code-block:: python
+      [('mass', 'none', 1.0f), ('length', 'none', 1.0f), ('time', 'none', 1.0f), ('G', 'none', 1.0f)]
+
+
+Units in pyEXP
+--------------
+
+The ``pyEXP.basis`` classes will use the units specified by the user
+in the ``pyEXP.coefs`` class created from simulation snapshots (see
+:ref:`intro-pyEXP-tutorial`) or read by the coefficient factory from
+an existing HDF5 file to produce accelerations using the value of the
+gravitational constant and length units provided by the user.
+
+
+.. _units-interface:
+
+The Units interface
+-------------------
+
+The `pyEXP` user interface includes two member functions for explicitly
+setting and removing units as part of the `Coef` class.  For setting
+units, we have:
+
+   .. code-block:: python
+
+      setUnits(type, unit, value)
+      setUnits(list)
+      removeUnits(type)
+      getAllowedUnitTypes()
+      getAllowedTypeAliases(type)
+      getAllowedUnitName(type)
+
+where ``type`` and ``unit`` are strings and ``value`` is a float.  The
+list is a list of tuples of ``(name, unit, value)``.  The last three
+members return the list of unit types, the recognized aliases for each
+type, and the allowed unit names for each unit type, respectively.
+
+For an example, suppose you are making a set of coefficients for a
+simulation with default Gadget units.  Say your coefficients instance
+is called ``coefs``.   The following command will register the unit set:
+
+   .. code-block:: python
+
+      coefs.setUnits([ ('mass', 'Msun', 1e10), ('length', 'kpc', 1.0),
+      ('velocity', 'km/s', 1.0), ('G', 'mixed', 43007.1) ])
+
+These units will be in the HDF5 that you create using
+:python:`coefs.WriteH5Coefs('filename')`.  You can query, for example,
+the allowed 'mass' units with the call
+:python:`coefs.getAllowedUnitNames('mass')` which returns:
+
+    .. code-block:: python
+
+       ['gram', 'earth_mass', 'solar_mass', 'kilograms', 'kg', 'g', 'Mearth', 'Msun', 'None', 'none']
+
+The allowed mass units for 'G' are:
+
+    .. code-block:: python
+
+       ['gram', 'earth_mass', 'solar_mass', 'kilograms', 'kg', 'g', 'Mearth', 'Msun', 'None', 'none']
+
+
+A quick note: 'mixed' is an allowed alias when dealing with
+gravitational constants that have physical units.  You can see all
+unit types with :python:`getAllowedUnitTypes()`; this returns
+:python:`['mass', 'length', 'time', 'velocity', 'G']`.  You can see
+the recognized aliases for each type using
+:python:`getAllowedTypeAliases(type)`.  The gravitational constant is
+	a special case. The recognized unit names for 'G' are: :python:`['','mixed', 'none', 'unitless']``.
+
+The C++ UI echos the functions above and adds functions to retrieve
+units
+
+   .. code-block:: C++
+
+      // Add or replace a unit by specifying its unit tuple
+      void setUnits(const std::string name, const std::string unit, float value);
+      // Add or replace a unit(s) by specifying an array of unit tuple
+      void setUnits(std::vector<std::tuple<std::string, std::string, float>> list);
+      // Remove a unit tuple by type
+      void removeUnits(const std::string type);
+      // Retrieve the currently define unit set
+      std::vector<std::tuple<std::string, std::string, float>> getUnits();
+      // Get a list of unit types and their aliases
+      std::vector<std::string> getAllowedTypes();
+      // Get a list aliases for each type
+      std::vector<std::string> getAllowedTypeAliases(const std::string& type);
+      // Get a list of unit name and their aliases for a given unit type
+      std::vector<std::string> getAllowedUnits(const std::string& type)
+
+and to interact with HDF files that will only be of interest to
+developers creating new coefficient classes.
+
+
+.. _hdf5-unit-attributes:
+
+The HDF5 specification
+----------------------
+
+.. note:: This and the following sections assumes basic familiarity
+          with HDF5 and `h5py` in particular.
+
+EXP HDF5 coefficient files contain a full set of metadata describing
+their basis type, basis parameters, geometry, and units as attributes
+of the root ("/") group.  The ``snapshots`` group contains all of the
+coefficient data and metadata (such as center and rotation
+information) for each snapshot time in the coefficient set.
+
+The units information is stored in the root group as dataset named
+"Units".  The dataset is a sequence  or list of 4 tuples.  Each tuple
+has three fields: two fixed length strings of sixteen (16) characters
+and a float value.
+
+For an EXP run, the units specification appears as dataset in the root
+group with the following hierarchy::
+
+     DATASET "Units" {
+      DATATYPE  H5T_COMPOUND {
+         H5T_STRING {
+            STRSIZE 16;
+            STRPAD H5T_STR_NULLTERM;
+            CSET H5T_CSET_UTF8;
+            CTYPE H5T_C_S1;
+         } "name";
+         H5T_STRING {
+            STRSIZE 16;
+            STRPAD H5T_STR_NULLTERM;
+            CSET H5T_CSET_UTF8;
+            CTYPE H5T_C_S1;
+         } "unit";
+         H5T_IEEE_F32LE "value";
+      }
+      DATASPACE  SIMPLE { ( 4 ) / ( 4 ) }
+      DATA {
+      (0): {
+            "G",
+            "none",
+            1
+         },
+      (1): {
+            "length",
+            "none",
+            1
+         },
+      (2): {
+            "mass",
+            "none",
+            1
+         },
+      (3): {
+            "time",
+            "none",
+            1
+         }
+      }
+   }
+
+
+.. _units_old:
+
+Backward compatibility
+----------------------
+
+A units attribute list could be easily added to existing HDF5 files
+using `h5py` using the schema described above.  For example, assuming
+that you already have an open HDF5 coefficient file named `f`, the
+units described in the scheme above could be added to `f` with the
+following code:
+
+   .. code-block:: python
+
+      import h5py
+      import numpy as np
+
+      # Define the compound datatype with fixed-length ASCII strings and a float32
+      dt = np.dtype([
+          ('name', 'S16'),   # Fixed-length ASCII string of 16 bytes
+          ('unit', 'S16'),   # Fixed-length ASCII string of 16 bytes
+          ('value', '<f4')   # Little-endian 32-bit float
+      ])
+
+      # Create the dataset data as a numpy array
+      data = np.array([
+          (b'G',       b'none', 1.0),
+          (b'length',  b'none', 1.0),
+          (b'mass',    b'none', 1.0),
+          (b'time',    b'none', 1.0),
+      ], dtype=dt)
+
+      # Open the HDF5 coefficient file and add the dataset
+      with h5py.File('outcoef.disk.runtag', 'a') as f:
+          # 'a' mode opens file for read/write, creates if not exist
+
+          # Create dataset 'Units'
+          if 'Units' in f:
+              # Optional: remove existing dataset if you want to overwrite
+              del f['Units']  
+          dset = f.create_dataset('Units', data=data, dtype=dt)
+
+      # File is automatically closed on leaving the 'with' block, flush and
+      # update saved data on disk
+
+You can update the  :python:`data` array in the example above to match your
+unit set.