WIP: Transforms

.gitignore

@@ -11,3 +11,4 @@ hexrd.egg-info
 *.*~
 *#
 .#*
+.pytest_cache

devel-conda-environment.yml

@@ -0,0 +1,22 @@
+name: hexrd-develop
+
+dependencies:
+  - python        =2
+  - dill
+  - fabio
+  - h5py
+  - numba
+  - numpy
+  - progressbar   >=2.3
+  - pyyaml
+  - scikit-learn
+  - scipy
+  - pytest
+  - ipython
+
+# UI related only?
+#  - matplotlib
+#  - python.app # [osx]
+#  - scikit-image
+#  - qtconsole
+#  - wxpython

hexrd/constants.py

@@ -66,6 +66,9 @@
 zeros_6x1 = np.zeros((6, 1))
 
 # reference beam direction and eta=0 ref in LAB FRAME for standard geometry
+# WARNING: In some parts of the code it is assumed that the beam rotation
+#          matrix is assumed to be the identity when using the standard
+#          beam_vec and eta_vec. If this changes code must be reviewed!
 beam_vec = -lab_z
 eta_vec = lab_x
 

hexrd/transforms/__init__.py

@@ -0,0 +1,33 @@
+"""Transforms module.
+
+Contains different implementations based on Python+Numpy, numba and
+a supporting C module. All three should adhere to the same interface,
+but performance will vary.
+
+Use the functions under this module scope to use the preferred versions,
+import the specific submodule if you want to use a specific version.
+"""
+from __future__ import absolute_import
+
+# While this may not be ideal, use an import line per API function.
+
+from .xf_numpy import angles_to_gvec
+from .xf_numpy import angles_to_dvec
+from .xf_numpy import gvec_to_xy
+from .xf_numpy import xy_to_gvec
+from .xf_numpy import solve_omega
+
+# utility functions
+from .xf_numpy import angular_difference
+from .xf_numpy import map_angle
+from .xf_numpy import row_norm
+from .xf_numpy import unit_vector
+from .xf_numpy import make_sample_rmat
+from .xf_numpy import make_rmat_of_expmap
+from .xf_numpy import make_binary_rmat
+from .xf_numpy import make_beam_rmat
+from .xf_numpy import angles_in_range
+from .xf_numpy import validate_angle_ranges
+from .xf_numpy import rotate_vecs_about_axis
+from .xf_numpy import quat_product_matrix
+from .xf_numpy import quat_distance

hexrd/transforms/transforms_CAPI.c → hexrd/transforms/src/transforms_CAPI.c

@@ -884,10 +884,11 @@ static PyObject * makeBinaryRotMat(PyObject * self, PyObject * args)
 
 static PyObject * makeEtaFrameRotMat(PyObject * self, PyObject * args)
 {
-  PyArrayObject *bHat, *eHat, *rMat;
+  PyArrayObject *bHat, *eHat, *rMat=NULL;
   int db, de;
   npy_intp nb, ne, dims[2];
   double *bPtr, *ePtr, *rPtr;
+  int errcode;
 
   /* Parse arguments */
   if ( !PyArg_ParseTuple(args,"OO", &bHat,&eHat)) return(NULL);
@@ -907,15 +908,35 @@ static PyObject * makeEtaFrameRotMat(PyObject * self, PyObject * args)
   dims[0] = 3; dims[1] = 3;
   rMat = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0);
 
+  if (rMat == NULL)
+      goto fail;
+
   /* Grab pointers to the various data arrays */
   bPtr = (double*)PyArray_DATA(bHat);
   ePtr = (double*)PyArray_DATA(eHat);
   rPtr = (double*)PyArray_DATA(rMat);
 
   /* Call the actual function */
-  makeEtaFrameRotMat_cfunc(bPtr,ePtr,rPtr);
+  errcode = makeEtaFrameRotMat_cfunc(bPtr, ePtr, rPtr);
+
+  if (errcode == 0) {
+      /* No error, return the matrix */
+      return((PyObject*)rMat);
+  } /* else ... fall back to fail code, but generating the appropriate error before */
+
+  switch (errcode) {
+  case TF_MAKE_BEAM_RMAT_ERR_BEAM_ZERO: /* beam vec is zero (within an epsilon) */
+      PyErr_SetString(PyExc_RuntimeError, "bvec_l MUST NOT be ZERO!");
+      break;
+  case TF_MAKE_BEAM_RMAT_ERR_COLLINEAR: /* beam vec and eta vec are collinear */
+      PyErr_SetString(PyExc_RuntimeError, "bvec_l and evec_l MUST NOT be collinear!");
+      break;
+  }
 
-  return((PyObject*)rMat);
+ fail:
+  /* Free the array if allocated, since it won't be returned */
+  Py_XDECREF(rMat);
+  return NULL;
 }
 
 static PyObject * validateAngleRanges(PyObject * self, PyObject * args)

hexrd/transforms/transforms_CFUNC.c → hexrd/transforms/src/transforms_CFUNC.c

@@ -11,6 +11,18 @@ static double sqrt_epsf = 1.5e-8;
 
 static double Zl[3] = {0.0,0.0,1.0};
 
+static inline
+void set_identity_matrix33(double * rPtr)
+{
+  int i;
+  for (i = 0; i < 9; i ++) {
+    if ( i%4 != 0 )
+      rPtr[i] = 0.0;
+    else
+      rPtr[i] = 1.0;
+  }
+}
+
 /******************************************************************************/
 /* Functions */
 
@@ -550,7 +562,10 @@ void oscillAnglesOfHKLs_cfunc(long int npts, double * hkls, double chi,
 /******************************************************************************/
 /* Utility Funtions */
 
-void unitRowVector_cfunc(int n, double * cIn, double * cOut)
+/* Stores normalized cIn into Cout. Cin and Cout being vectors of n elements
+ * returns 0 if normalization was possible, 1 if the normalization wasn't
+ */
+int unitRowVector_cfunc(int n, double * cIn, double * cOut)
 {
   int j;
   double nrm;
@@ -564,10 +579,12 @@ void unitRowVector_cfunc(int n, double * cIn, double * cOut)
     for (j=0; j<n; j++) {
       cOut[j] = cIn[j]/nrm;
     }
+    return 0;
   } else {
     for (j=0; j<n; j++) {
       cOut[j] = cIn[j];
     }
+    return 1;
   }
 }
 
@@ -637,15 +654,9 @@ void makeOscillRotMat_cfunc(double chi, double ome, double * rPtr)
 
 void makeRotMatOfExpMap_cfunc(double * ePtr, double * rPtr)
 {
-  int i;
   double c, s, phi;
 
-  for (i=0; i<9; i++) {
-    if ( i%4 != 0 )
-      rPtr[i] = 0.0;
-    else
-      rPtr[i] = 1.0;
-  }
+  set_identity_matrix33(rPtr);
 
   phi = sqrt(ePtr[0]*ePtr[0]+ePtr[1]*ePtr[1]+ePtr[2]*ePtr[2]);
 
@@ -675,7 +686,7 @@ void makeRotMatOfExpMap_cfunc(double * ePtr, double * rPtr)
 
 void makeRotMatOfQuat_cfunc(int nq, double * qPtr, double * rPtr)
 {
-  int i, j;
+  int i;
   double c, s, phi, n[3]={0.0,0.0,0.0};
 
   for (i=0; i<nq; i++) {
@@ -700,12 +711,7 @@ void makeRotMatOfQuat_cfunc(int nq, double * qPtr, double * rPtr)
       rPtr[9*i+8] = c + n[2]*n[2]*(1. - c);
     }
     else {
-      for (j=0; j<9; i++) {
-	if ( j%4 == 0 )
-	  rPtr[9*i+j] = 1.0;
-	else
-	  rPtr[9*i+j] = 0.0;
-      }
+      set_identity_matrix33(rPtr+9*i);
     }
   }
 }
@@ -722,25 +728,40 @@ void makeBinaryRotMat_cfunc(double * aPtr, double * rPtr)
   }
 }
 
-void makeEtaFrameRotMat_cfunc(double * bPtr, double * ePtr, double * rPtr)
+int 
+makeEtaFrameRotMat_cfunc(double * bPtr, double * ePtr, double * rPtr)
 {
   /*
    * This function generates a COB matrix that takes components in the
    * BEAM frame to LAB frame
    *
    * NOTE: the beam and eta vectors MUST NOT BE COLINEAR!!!!
+   * 
+   * Added a return value acting as an error code. Possible values are:
+   * 0 - No error
+   * 1 - beam vector is zero within an epsilon
+   * 2 - beam and eta vectors are collinear.
    */
-  int i;
+  int i, err;
   double yPtr[3], bHat[3], yHat[3], xHat[3];
 
+  err = unitRowVector_cfunc(3, bPtr, bHat);
+  if (0 != err) {
+      /* can't normalize beam vector due to being zero */
+      return TF_MAKE_BEAM_RMAT_ERR_BEAM_ZERO;
+  }
+
   /* find Y as e ^ b */
   yPtr[0] = ePtr[1]*bPtr[2] - bPtr[1]*ePtr[2];
   yPtr[1] = ePtr[2]*bPtr[0] - bPtr[2]*ePtr[0];
   yPtr[2] = ePtr[0]*bPtr[1] - bPtr[0]*ePtr[1];
 
-  /* Normalize beam (-Z) and Y vectors */
-  unitRowVector_cfunc(3,bPtr,bHat);
-  unitRowVector_cfunc(3,yPtr,yHat);
+  /* Normalize e ^ b */
+  err = unitRowVector_cfunc(3,yPtr,yHat);
+  if (0 != err) {
+      /* e ^ b is close to zero... eta and beam are collinear... */
+      return TF_MAKE_BEAM_RMAT_ERR_COLLINEAR;
+  }
 
   /* Find X as b ^ Y */
   xHat[0] = bHat[1]*yHat[2] - yHat[1]*bHat[2];
@@ -753,6 +774,8 @@ void makeEtaFrameRotMat_cfunc(double * bPtr, double * ePtr, double * rPtr)
     rPtr[3*i+1] = yHat[i];
     rPtr[3*i+2] = -bHat[i];
   }
+
+  return 0; /* no error */
 }
 
 void validateAngleRanges_old_cfunc(int na, double * aPtr, int nr, double * minPtr, double * maxPtr, bool * rPtr)

hexrd/transforms/transforms_CFUNC.h → hexrd/transforms/src/transforms_CFUNC.h

@@ -61,7 +61,8 @@ void oscillAnglesOfHKLs_cfunc(long int npts, double * hkls, double chi,
 /******************************************************************************/
 /* Utility Funtions */
 
-void unitRowVector_cfunc(int n, double * cIn, double * cOut);
+int 
+unitRowVector_cfunc(int n, double * cIn, double * cOut);
 
 void unitRowVectors_cfunc(int m, int n, double * cIn, double * cOut);
 
@@ -75,7 +76,11 @@ void makeRotMatOfQuat_cfunc(int nq, double * qPtr, double * rPtr);
 
 void makeBinaryRotMat_cfunc(double * aPtr, double * rPtr);
 
-void makeEtaFrameRotMat_cfunc(double * bPtr, double * ePtr, double * rPtr);
+
+#define TF_MAKE_BEAM_RMAT_ERR_BEAM_ZERO 1
+#define TF_MAKE_BEAM_RMAT_ERR_COLLINEAR 2
+int 
+makeEtaFrameRotMat_cfunc(double * bPtr, double * ePtr, double * rPtr);
 
 void validateAngleRanges_cfunc(int na, double * aPtr, int nr, double * minPtr, double * maxPtr, bool * rPtr, int ccw);
 

hexrd/transforms/tests/README.md

@@ -0,0 +1,91 @@
+# Tests for the transforms module #
+
+This directory contains the tests for the transforms module. You can
+run them by running `pytest <this directory>` in the command line.
+
+
+## Physical structure of the transform package unit tests ##
+
+The tests are organized in a test file for each function supported in
+the API. The test file is named `test_<function>.py`.
+
+The tests are parametrized so that the same tests are run in all
+available implementations. This is achieved by using
+`pytest.mark.parametrize`.  Typically, a decorator `all_impls` (short
+for "all implementations") is defined as a `pytest.mark.parametrize`
+that define two arguments, the function implementation and a submodule
+name (for reporting mostly). A list of parametrizations defines the
+implementations to test. The decorator will be added to each test
+function so that the test is run on all the different implementations.
+
+An example, for an API function named `foobar`, having a numpy and a numba
+implementation (but not capi) would be:
+
+```python
+from .. import foobar as default_foobar
+from ..xf_numpy import foobar as numpy_foobar
+from ..xf_numba import foobar as numba_foobar
+
+all_impls = pytest.mark.parametrize('foobar_impl', 'module_name',
+	[(numpy_foobar, 'numpy'),
+	 (numba_foobar, 'numba'),
+	 (default_foobar, 'default')]
+)
+
+@allimpls
+def test_sample_test(foobar_impl, module_name):
+	<Your test goes here>
+```
+
+In the example, the test `test_sample_test` will be run thrice. First
+using the implementation of `foobar` in `xf_numpy`. Then using the
+implementation of `foobar` in `xf_numba`. Finally using the
+implementation exported at the module level.
+
+Typically the implementation exported at the module level will be one
+of the implementations explicitly tested from one of the
+implementation packages. In that case the tests on that implementation
+will be run twice, once as the module level export and again for the
+submodule implementation.
+
+
+## Recommented tests ##
+
+When writing the unit tests for transforms functions, testing the
+following features should be considered:
+
+* Test basic functionality, preferably using trivial/known setups as
+  arguments.
+
+* Test default arguments against their equivalent explicit versions.
+  In many cases the default argument for array arguments is `None`,
+  but semantically a `None` value is logically equivalent to some
+  array constant. The explicit check should use that constant, not
+  `None`.
+
+* Test error cases that should raise exceptions, if they exist. Check
+  that the exception happens. If the numerical code has some kind of
+  singularity, even if it does not raise exceptions, write a test that
+  checks the singularity is handled in a consistent (and preferably,
+  documented) way.
+
+* If the function broadcasts, test that it broadcasts properly. But
+  check also the arguments that do not broadcast. Some implementations
+  may use different code paths for broadcasting and scalar versions.
+
+* If a function broadcasts and it has some kind of singularity, make
+  sure that tests exist checking arguments where a singularity happens
+  while broadcasting, mixed with cases where the singularity doesn't
+  happen.
+
+* Functions taking numpy arrays as arguments should be tested with
+  arguments that are strided in a non natural way. Implementations
+  based on C code are prone to not handling those properly if not
+  carefully implemented/wrapped.
+
+Testing behavior on handling special values on input, like `NaN`, may
+be important in some functions.
+
+Finally, when a bug is found, it is desirable to add a test that
+reproduces that bug and leave it so that any regression does not go
+unnoticed.

hexrd/transforms/tests/__init__.py

@@ -0,0 +1,17 @@
+"""Test suite for hexrd.transforms"""
+
+# Use a file for each api function to test named test_<function>.py
+#
+# Configure the all_impls decorator (based on pytest.mark.parametrize)
+#
+# Each test to be applied to all implementations should be preceded by an
+# @all_impls decorator. The test function will take two arguments, the first one
+# is the implementation to test, the second will be a "module" name used to
+# better identify the function implementation.
+#
+# It should be possible to test different implementations in the same module
+# that differ in name, and distinguish between implementation from different
+# modules that use the same name by appropriate tagging with a module name.
+
+
+