From 017c6c56d155a90e63558492be770155ce7a7f58 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 23 Apr 2018 15:35:57 +0200 Subject: [PATCH 01/30] added a conda development environment yml file --- devel-conda-environment.yml | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 devel-conda-environment.yml diff --git a/devel-conda-environment.yml b/devel-conda-environment.yml new file mode 100644 index 00000000..4cde6896 --- /dev/null +++ b/devel-conda-environment.yml @@ -0,0 +1,21 @@ +name: hexrd-develop + +dependencies: + - python =2 + - dill + - fabio + - h5py + - numba + - numpy + - progressbar >=2.3 + - pyyaml + - scikit-learn + - scipy + - ipython + +# UI related only? +# - matplotlib +# - python.app # [osx] +# - scikit-image +# - qtconsole +# - wxpython From c723e513c797fb47bf1255dc14316050fdc6c417 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 23 Apr 2018 15:52:03 +0200 Subject: [PATCH 02/30] fixed possible infinite loop on makeRotMatOfQuat_cfunc --- hexrd/transforms/transforms_CFUNC.c | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/hexrd/transforms/transforms_CFUNC.c b/hexrd/transforms/transforms_CFUNC.c index 8cb9f7d8..14add20f 100644 --- a/hexrd/transforms/transforms_CFUNC.c +++ b/hexrd/transforms/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 */ @@ -640,12 +652,7 @@ 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]); @@ -700,12 +707,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); } } } From d6049c4af88abfb6b57cfa1c5d36bf4cf0c8e86a Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 23 Apr 2018 16:22:57 +0200 Subject: [PATCH 03/30] removed a couple of variable unused warnings --- hexrd/transforms/transforms_CFUNC.c | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/hexrd/transforms/transforms_CFUNC.c b/hexrd/transforms/transforms_CFUNC.c index 14add20f..7f15f355 100644 --- a/hexrd/transforms/transforms_CFUNC.c +++ b/hexrd/transforms/transforms_CFUNC.c @@ -649,7 +649,6 @@ void makeOscillRotMat_cfunc(double chi, double ome, double * rPtr) void makeRotMatOfExpMap_cfunc(double * ePtr, double * rPtr) { - int i; double c, s, phi; set_identity_matrix33(rPtr); @@ -682,7 +681,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 Date: Mon, 23 Apr 2018 16:45:33 +0200 Subject: [PATCH 04/30] added a minimal __init__.py for the transforms module. --- hexrd/transforms/__init__.py | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 hexrd/transforms/__init__.py diff --git a/hexrd/transforms/__init__.py b/hexrd/transforms/__init__.py new file mode 100644 index 00000000..641f9c12 --- /dev/null +++ b/hexrd/transforms/__init__.py @@ -0,0 +1,11 @@ +"""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. +""" + +# TODO: make public the default definitions From 3738d7222a9ef177c8ed26a1f5c1065afd141274 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 23 Apr 2018 17:33:59 +0200 Subject: [PATCH 05/30] moved transforms C code into src subdirectory. Modified setup.py accordingly. Also made coherent the location of the sglite library location with that of transforms. --- hexrd/transforms/{ => src}/Makefile | 0 hexrd/transforms/{ => src}/stdbool.h | 0 hexrd/transforms/{ => src}/transforms_CAPI.c | 0 hexrd/transforms/{ => src}/transforms_CAPI.h | 0 hexrd/transforms/{ => src}/transforms_CFUNC.c | 0 hexrd/transforms/{ => src}/transforms_CFUNC.h | 0 hexrd/xrd/spacegroup.py | 2 +- hexrd/xrd/transforms_CAPI.py | 2 +- setup.py | 16 ++++++++-------- 9 files changed, 10 insertions(+), 10 deletions(-) rename hexrd/transforms/{ => src}/Makefile (100%) rename hexrd/transforms/{ => src}/stdbool.h (100%) rename hexrd/transforms/{ => src}/transforms_CAPI.c (100%) rename hexrd/transforms/{ => src}/transforms_CAPI.h (100%) rename hexrd/transforms/{ => src}/transforms_CFUNC.c (100%) rename hexrd/transforms/{ => src}/transforms_CFUNC.h (100%) diff --git a/hexrd/transforms/Makefile b/hexrd/transforms/src/Makefile similarity index 100% rename from hexrd/transforms/Makefile rename to hexrd/transforms/src/Makefile diff --git a/hexrd/transforms/stdbool.h b/hexrd/transforms/src/stdbool.h similarity index 100% rename from hexrd/transforms/stdbool.h rename to hexrd/transforms/src/stdbool.h diff --git a/hexrd/transforms/transforms_CAPI.c b/hexrd/transforms/src/transforms_CAPI.c similarity index 100% rename from hexrd/transforms/transforms_CAPI.c rename to hexrd/transforms/src/transforms_CAPI.c diff --git a/hexrd/transforms/transforms_CAPI.h b/hexrd/transforms/src/transforms_CAPI.h similarity index 100% rename from hexrd/transforms/transforms_CAPI.h rename to hexrd/transforms/src/transforms_CAPI.h diff --git a/hexrd/transforms/transforms_CFUNC.c b/hexrd/transforms/src/transforms_CFUNC.c similarity index 100% rename from hexrd/transforms/transforms_CFUNC.c rename to hexrd/transforms/src/transforms_CFUNC.c diff --git a/hexrd/transforms/transforms_CFUNC.h b/hexrd/transforms/src/transforms_CFUNC.h similarity index 100% rename from hexrd/transforms/transforms_CFUNC.h rename to hexrd/transforms/src/transforms_CFUNC.h diff --git a/hexrd/xrd/spacegroup.py b/hexrd/xrd/spacegroup.py index 1cad811a..4bda6d6d 100644 --- a/hexrd/xrd/spacegroup.py +++ b/hexrd/xrd/spacegroup.py @@ -73,7 +73,7 @@ from collections import OrderedDict from math import sqrt, floor -import hexrd.xrd.sglite as sglite +import hexrd.sglite as sglite # __all__ = ['SpaceGroup'] # diff --git a/hexrd/xrd/transforms_CAPI.py b/hexrd/xrd/transforms_CAPI.py index 1eb71f13..6eb1e4d7 100644 --- a/hexrd/xrd/transforms_CAPI.py +++ b/hexrd/xrd/transforms_CAPI.py @@ -29,7 +29,7 @@ import numpy as np import sys -from hexrd.xrd import _transforms_CAPI +from hexrd import _transforms_CAPI from numpy import float_ as nFloat from numpy import int_ as nInt diff --git a/setup.py b/setup.py index cb8bfd5c..5c0190f2 100644 --- a/setup.py +++ b/setup.py @@ -66,26 +66,26 @@ def run(self): # for SgLite -srclist = [ +sglite_srcs = [ 'sgglobal.c', 'sgcb.c', 'sgcharmx.c', 'sgfile.c', 'sggen.c', 'sghall.c', 'sghkl.c', 'sgltr.c', 'sgmath.c', 'sgmetric.c', 'sgnorm.c', 'sgprop.c', 'sgss.c', 'sgstr.c', 'sgsymbols.c', 'sgtidy.c', 'sgtype.c', 'sgutil.c', 'runtests.c', 'sglitemodule.c' ] -srclist = [os.path.join('hexrd/sglite', f) for f in srclist] +sglite_srcs = [os.path.join('hexrd','sglite', f) for f in sglite_srcs] sglite_mod = Extension( - 'hexrd.xrd.sglite', - sources=srclist, + 'hexrd.sglite', + sources=sglite_srcs, define_macros=[('PythonTypes', 1)] ) # for transforms -srclist = ['transforms_CAPI.c', 'transforms_CFUNC.c'] -srclist = [os.path.join('hexrd/transforms', f) for f in srclist] +transforms_srcs = ['transforms_CAPI.c', 'transforms_CFUNC.c'] +transforms_srcs = [os.path.join('hexrd', 'transforms', 'src' , f) for f in transforms_srcs] transforms_mod = Extension( - 'hexrd.xrd._transforms_CAPI', - sources=srclist, + 'hexrd._transforms_CAPI', + sources=transforms_srcs, include_dirs=[np_include_dir] ) From 45414c5ca729394081e1f14ea3dfb498436f43d8 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 23 Apr 2018 17:34:29 +0200 Subject: [PATCH 06/30] added pytest to the devel-conda-environment yaml --- devel-conda-environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/devel-conda-environment.yml b/devel-conda-environment.yml index 4cde6896..2ab7b85d 100644 --- a/devel-conda-environment.yml +++ b/devel-conda-environment.yml @@ -11,6 +11,7 @@ dependencies: - pyyaml - scikit-learn - scipy + - pytest - ipython # UI related only? From 90fc80bf687dcbee146c2747bbb25730dfbdfbda Mon Sep 17 00:00:00 2001 From: Donald Boyce Date: Mon, 7 May 2018 11:41:52 -0400 Subject: [PATCH 07/30] name and path changes to get xf_numpy to import --- hexrd/transforms/xf_numpy.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/hexrd/transforms/xf_numpy.py b/hexrd/transforms/xf_numpy.py index c34668f6..9cdb549e 100644 --- a/hexrd/transforms/xf_numpy.py +++ b/hexrd/transforms/xf_numpy.py @@ -32,7 +32,7 @@ from numpy import int_ as npint # from hexrd import constants as cnst -import constants as cnst +import hexrd.constants as cnst # ============================================================================= @@ -43,7 +43,7 @@ # and multiply through? def _beam_to_crystal(vecs, rmat_b=None, rmat_s=None, rmat_c=None): """ - Helper function to take vectors definced in the BEAM frame through LAB + Helper function to take vectors defined in the BEAM frame through LAB to either SAMPLE or CRYSTAL """ @@ -208,7 +208,7 @@ def _z_project(x, y): def angles_to_gvec( angs, - beam_vec=cnst.ref_beam_vec, eta_vec=cnst.ref_eta_vec, + beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, chi=None, rmat_c=None): """ Takes triplets of angles in the beam frame (2*theta, eta, omega) @@ -235,7 +235,7 @@ def angles_to_gvec( def angles_to_dvec( angs, - beam_vec=cnst.ref_beam_vec, eta_vec=cnst.ref_eta_vec, + beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, chi=None, rmat_c=None): """ Takes triplets of angles in the beam frame (2*theta, eta, omega) @@ -263,7 +263,7 @@ def angles_to_dvec( def gvec_to_xy(gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, - beam_vec=cnst.ref_beam_vec, + beam_vec=cnst.beam_vec, vmat_inv=None, bmat=None): """ @@ -456,7 +456,7 @@ def xy_to_gvec(xy_d, npts = len(xy_d) # need beam vector - bhat_l = cnst.ref_beam_vec + bhat_l = cnst.beam_vec if rmat_b is not None: bhat_l = -rmat_b[:, 2] @@ -513,7 +513,7 @@ def solve_omega(gvecs, chi, rmat_c, wavelength, (3, 3) COB matrix taking components in the CRYSTAL FRAME to the SAMPLE FRAME wavelength : float - The X-ray wavelength in Ångstroms + The X-ray wavelength in Angstroms bmat : array_like, optional The (3, 3) COB matrix taking components in the RECIPROCAL LATTICE FRAME to the CRYSTAL FRAME; if supplied, it is From f8d95851b34911c477331c8cadfb13b07d59116e Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 14 May 2018 12:13:11 +0200 Subject: [PATCH 08/30] api definitions and decorator, as well as a minimal __init__ --- hexrd/transforms/__init__.py | 24 +- hexrd/transforms/transforms_definitions.py | 540 +++++++++++++++++++++ hexrd/transforms/xf_numpy.py | 334 +------------ 3 files changed, 586 insertions(+), 312 deletions(-) create mode 100644 hexrd/transforms/transforms_definitions.py diff --git a/hexrd/transforms/__init__.py b/hexrd/transforms/__init__.py index 641f9c12..65ae3701 100644 --- a/hexrd/transforms/__init__.py +++ b/hexrd/transforms/__init__.py @@ -7,5 +7,27 @@ 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 -# TODO: make public the default definitions +# 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 diff --git a/hexrd/transforms/transforms_definitions.py b/hexrd/transforms/transforms_definitions.py new file mode 100644 index 00000000..988eaf6f --- /dev/null +++ b/hexrd/transforms/transforms_definitions.py @@ -0,0 +1,540 @@ +"""This module provides the definitions for the transforms API. It will also +provide a decorator to add to any implementation of the API. This module will +contain the reference comment that will be added to any function that implements +an API function, as well as a means to add pre and post conditions as an additional +way to document the implementations. + +Pre and Post conditions will be in the form of code, there will be means to +execute the scripts forcing those conditions to be evaluated and raise errors +if they are not met. This should always be optional and incur on no overhead +unless enabled, only to be used for debugging and validation purposes. +""" + +import os +import functools + +CHECK_API = os.getenv("HEXRD_XF_CHECK") + +class DEF_Func(object): + """Documentation to use for the function""" + + @property + def _args(self): + """The position arguments the function is supposed to have""" + return None + + @property + def _kwargs(self): + """The keyword arguments the function is supposed to have""" + return None + + def _PRECOND(*arg, **kwarg): + pass + + def _POSTCOND(results, *args, **kwargs): + pass + + +# ============================================================================== +# API +# ============================================================================== + +class DEF_angles_to_gvec(DEF_Func): + """ + Takes triplets of angles in the beam frame (2*theta, eta, omega) + to components of unit G-vectors in the LAB frame. If the omega + values are not trivial (i.e. angs[:, 2] = 0.), then the components + are in the SAMPLE frame. If the crystal rmat is specified and + is not the identity, then the components are in the CRYSTAL frame. + """ + @property + def _args(self): + return ('angs',) + + @property + def _kwargs(self): + return ('beam_vec', 'eta_vec', 'chi', 'rmat_c') + + +class DEF_angles_to_dvec(DEF_Func): + """ + Takes triplets of angles in the beam frame (2*theta, eta, omega) + to components of unit diffraction vectors in the LAB frame. If the + omega values are not trivial (i.e. angs[:, 2] = 0.), then the + components are in the SAMPLE frame. If the crystal rmat is specified + and is not the identity, then the components are in the CRYSTAL frame. + """ + @property + def _args(self): + return ('angs',) + + @property + def _kwargs(self): + return ('beam_vec', 'eta_vec', 'chi', 'rmat_c') + + +class DEF_gvec_to_xy(DEF_Func): + """ + Takes a concatenated list of reciprocal lattice vectors components in the + CRYSTAL FRAME to the specified detector-relative frame, subject to the + following: + + 1) it must be able to satisfy a bragg condition + 2) the associated diffracted beam must intersect the detector plane + + Parameters + ---------- + gvec_c : array_like + Concatenated triplets of G-vector components in the CRYSTAL FRAME. + rmat_d : array_like + The (3, 3) COB matrix taking components in the + DETECTOR FRAME to the LAB FRAME + rmat_s : array_like + The (3, 3) COB matrix taking components in the + SAMPLE FRAME to the LAB FRAME + rmat_c : array_like + The (3, 3) COB matrix taking components in the + CRYSTAL FRAME to the SAMPLE FRAME + tvec_d : array_like + The (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME + tvec_s : array_like + The (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME + tvec_c : array_like + The (3, ) translation vector connecting SAMPLE FRAME to CRYSTAL FRAME + beam_vec : array_like, optional + The (3, ) incident beam propagation vector components in the LAB FRAME; + the default is [0, 0, -1], which is the standard setting. + vmat_inv : array_like, optional + The (3, 3) matrix of inverse stretch tensor components in the + SAMPLE FRAME. The default is None, which implies a strain-free state + (i.e. V = I). + bmat : array_like, optional + The (3, 3) COB matrix taking components in the + RECIPROCAL LATTICE FRAME to the CRYSTAL FRAME; if supplied, it is + assumed that the input `gvecs` are G-vector components in the + RECIPROCL LATTICE FRAME (the default is None, which implies components + in the CRYSTAL FRAME) + + Returns + ------- + array_like + The (n, 2) array of [x, y] diffracted beam intersections for each of + the n input G-vectors in the DETECTOR FRAME (all Z_d coordinates are 0 + and excluded). For each input G-vector that cannot satisfy a Bragg + condition or intersect the detector plane, [NaN, Nan] is returned. + + Raises + ------ + AttributeError + The ``Raises`` section is a list of all exceptions + that are relevant to the interface. + ValueError + If `param2` is equal to `param1`. + + Notes + ----- + + """ + @property + def _args(self): + return ('gvec_c', 'rmat_d', 'rmat_s', 'rmat_c', 'tvec_d', 'tvec_s', + 'tvec_c') + + @property + def _kwargs(self): + return ('beam_vec', 'vmat_inv', 'bmat') + + +class DEF_xy_to_gvec(DEF_Func): + """ + Takes a list cartesian (x, y) pairs in the DETECTOR FRAME and calculates + the associated reciprocal lattice (G) vectors and (bragg angle, azimuth) + pairs with respect to the specified beam and azimth (eta) reference + directions. + + Parameters + ---------- + xy_d : array_like + (n, 2) array of n (x, y) coordinates in DETECTOR FRAME + rmat_d : array_like + (3, 3) COB matrix taking components in the + DETECTOR FRAME to the LAB FRAME + rmat_s : array_like + (3, 3) COB matrix taking components in the + SAMPLE FRAME to the LAB FRAME + tvec_d : array_like + (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME + tvec_s : array_like + (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME + tvec_c : array_like + (3, ) translation vector connecting SAMPLE FRAME to CRYSTAL FRAME + rmat_b : array_like, optional + (3, 3) COB matrix taking components in the BEAM FRAME to the LAB FRAME; + defaults to None, which implies the standard setting of identity. + distortion : distortion class, optional + Default is None + output_ref : bool, optional + If True, prepends the apparent bragg angle and azimuth with respect to + the SAMPLE FRAME (ignoring effect of non-zero tvec_c) + + Returns + ------- + array_like + (n, 2) ndarray containing the (tth, eta) pairs associated with each + (x, y) associated with gVecs + array_like + (n, 3) ndarray containing the associated G vector directions in the + LAB FRAME + array_like, optional + if output_ref is True + + Notes + ----- + ???: is there a need to flatten the tvec inputs? + ???: include optional wavelength input for returning G with magnitude? + ???: is there a need to check that rmat_b is orthogonal if spec'd? + """ + + @property + def _args(self): + return ('xy_d', 'rmat_d', 'rmat_s', 'tvec_d', 'tvec_s', 'tvec_c') + + @property + def _kwargs(self): + return ('rmat_b', 'distortion', 'output_ref') + + +class DEF_solve_omega(DEF_Func): + """ + For the monochromatic rotation method. + + Solve the for the rotation angle pairs that satisfy the bragg conditions + for an input list of G-vector components. + + Parameters + ---------- + gvecs : array_like + Concatenated triplets of G-vector components in either the + CRYSTAL FRAME or RECIPROCAL FRAME (see optional kwarg `bmat` below). + The shape when cast as a 2-d ndarray is (n, 3), representing n vectors. + chi : float + The inclination angle of the goniometer axis (standard coords) + rmat_c : array_like + (3, 3) COB matrix taking components in the + CRYSTAL FRAME to the SAMPLE FRAME + wavelength : float + The X-ray wavelength in Angstroms + bmat : array_like, optional + The (3, 3) COB matrix taking components in the + RECIPROCAL LATTICE FRAME to the CRYSTAL FRAME; if supplied, it is + assumed that the input `gvecs` are G-vector components in the + RECIPROCL LATTICE FRAME (the default is None, which implies components + in the CRYSTAL FRAME) + vmat_inv : array_like, optional + The (3, 3) matrix of inverse stretch tensor components in the + SAMPLE FRAME. The default is None, which implies a strain-free state + (i.e. V = I). + rmat_b : array_like, optional + (3, 3) COB matrix taking components in the BEAM FRAME to the LAB FRAME; + defaults to None, which implies the standard setting of identity. + + Returns + ------- + ome0 : array_like + The (n, 3) ndarray containing the feasible (tth, eta, ome) triplets for + each input hkl (first solution) + ome1 : array_like + The (n, 3) ndarray containing the feasible (tth, eta, ome) triplets for + each input hkl (second solution) + + Notes + ----- + The reciprocal lattice vector, G, will satisfy the the Bragg condition + when: + + b.T * G / ||G|| = -sin(theta) + + where b is the incident beam direction (k_i) and theta is the Bragg + angle consistent with G and the specified wavelength. The components of + G in the lab frame in this case are obtained using the crystal + orientation, Rc, and the single-parameter oscillation matrix, Rs(ome): + + Rs(ome) * Rc * G / ||G|| + + The equation above can be rearranged to yeild an expression of the form: + + a*sin(ome) + b*cos(ome) = c + + which is solved using the relation: + + a*sin(x) + b*cos(x) = sqrt(a**2 + b**2) * sin(x + alpha) + + --> sin(x + alpha) = c / sqrt(a**2 + b**2) + + where: + + alpha = arctan2(b, a) + + The solutions are: + + / + | arcsin(c / sqrt(a**2 + b**2)) - alpha + x = < + | pi - arcsin(c / sqrt(a**2 + b**2)) - alpha + \ + + There is a double root in the case the reflection is tangent to the + Debye-Scherrer cone (c**2 = a**2 + b**2), and no solution if the + Laue condition cannot be satisfied (filled with NaNs in the results + array here) + """ + @property + def _args(self): + return ('gvecs', 'chi', 'rmat_c', 'wavelength') + + @property + def _kwargs(self): + return ('bmat', 'vmat_inv', 'rmat_b') + + +# ============================================================================== +# UTILITY FUNCTIONS API +# ============================================================================== + +class DEF_angular_difference(DEF_Func): + """ + Do the proper (acute) angular difference in the context of a branch cut. + + *) Default angular range is [-pi, pi] + """ + @property + def _args(self): + return ('ang_list0', 'ang_list1') + + @property + def _kwargs(self): + return ('units',) + + +class DEF_map_angle(DEF_Func): + """ + Utility routine to map an angle into a specified period + + actual function is map_angle(ang[, range], units=cnst.angular_units). + range is optional and defaults to the appropriate angle for the unit + centered on 0. + """ + @property + def _args(self): + return ('ang',) + + @property + def _kwargs(self): + return ('range', 'units') + + +class DEF_row_norm(DEF_Func): + """ + normalize array of row vectors (vstacked, axis = 1) + """ + @property + def _args(self): + return ('a') + + +class DEF_unit_vector(DEF_Func): + """ + normalize an array of row vectors (vstacked, axis=0) + """ + @property + def _args(self): + return('a') + + +class DEF_make_sample_rmat(DEF_Func): + """ + Make SAMPLE frame rotation matrices as composition of + rotation of ome about the axis + + [0., cos(chi), sin(chi)] + + in the lab frame + """ + @property + def _args(self): + return ('chi', 'ome') + + +class DEF_make_rmat_of_expmap(DEF_Func): + """ + Calculates the rotation matrix from an exponential map + """ + @property + def _args(self): + return ('exp_map',) + + +class DEF_make_binary_rmat(DEF_Func): + """ + make a binary rotation matrix about the specified axis + """ + @property + def _args(self): + return ('n',) + + +class DEF_make_beam_rmat(DEF_Func): + """ + make eta basis COB matrix with beam antiparallel with Z + + takes components from BEAM frame to LAB + """ + @property + def _args(self): + return ('bvec_l', 'evec_l') + + +class DEF_angles_in_range(DEF_Func): + """Determine whether angles lie in or out of specified ranges + + *angles* - a list/array of angles + *starts* - a list of range starts + *stops* - a list of range stops + + OPTIONAL ARGS: + *degrees* - [True] angles & ranges in degrees (or radians) + """ + @property + def _args(self): + return ('angles', 'starts', 'stops') + + @property + def _kwargs(self): + return ('degrees',) + + +class DEF_validate_angle_ranges(DEF_Func): + """ + A better way to go. find out if an angle is in the range + CCW or CW from start to stop + There is, of course, an ambigutiy if the start and stop angle are + the same; we treat them as implying 2*pi having been mapped + """ + @property + def _args(self): + return ('ang_list', 'startAngs', 'stopAngs') + + @property + def _kwargs(self): + return ('ccw',) + + +class DEF_rotate_vecs_about_axis(DEF_Func): + """ + Rotate vectors about an axis + + INPUTS + *angle* - array of angles (len == 1 or n) + *axis* - array of unit vectors (shape == (3, 1) or (3, n)) + *vec* - array of vectors to be rotated (shape = (3, n)) + + Quaternion formula: + if we split v into parallel and perpedicular components w.r.t. the + axis of quaternion q, + + v = a + n + + then the action of rotating the vector dot(R(q), v) becomes + + v_rot = (q0**2 - |q|**2)(a + n) + 2*dot(q, a)*q + 2*q0*cross(q, n) + + """ + @property + def _args(self): + return ('angle', 'axis', 'vecs') + + +class DEF_quat_product_matrix(DEF_Func): + """ + Form 4 x 4 array to perform the quaternion product + + USAGE + qmat = quatProductMatrix(q, mult='right') + + INPUTS + 1) quats is (4,), an iterable representing a unit quaternion + horizontally concatenated + 2) mult is a keyword arg, either 'left' or 'right', denoting + the sense of the multiplication: + + / quatProductMatrix(h, mult='right') * q + q * h --> < + \ quatProductMatrix(q, mult='left') * h + + OUTPUTS + 1) qmat is (4, 4), the left or right quaternion product + operator + + NOTES + *) This function is intended to replace a cross-product based + routine for products of quaternions with large arrays of + quaternions (e.g. applying symmetries to a large set of + orientations). + """ + @property + def _args(self): + return ('q',) + + @property + def _kwargs(self): + return ('mult',) + + +class DEF_quat_distance(DEF_Func): + """ + find the distance between two unit quaternions under symmetry group + """ + @property + def _args(self): + return ('q1', 'q2', 'qsym') + + +# ============================================================================== +# Decorator to mark implementations of the API. Names must match. +# ============================================================================== + +def xf_api(f): + """decorator to apply to the entry points of the transforms module""" + api_call = f.__name__ + + try: + fn_def = globals()['DEF_'+api_call] + except KeyError: + raise RuntimeError("xf_api function '%s' doesn't have a definition.") + + try: + if not (isinstance(fn_def.__doc__, basestring) and + callable(fn_def._PRECOND) and + callable(fn_def._POSTCOND)): + raise Exception() + except Exception: + raise RuntimeError("xf_api definition for function '%s' seems incorrect.") + + # At this point use a wrapper that calls pre and post conditions if checking + # is enabled, otherwise leave the function "as is". + if CHECK_API: + @functools.wraps(f, assigned={__doc__: fn_def.__doc__}) + def wrapper(*args, **kwargs): + fn_def._PRECOND(*args, **kwargs) + result = f(*args, **kwargs) + fn_def._POSTCOND(result, *args, **kwargs) + return result + + return wrapper + else: + # just try to put the right documentation on the function + try: + f.__doc__ = fn_def.__doc__ + except Exception: + pass + return f diff --git a/hexrd/transforms/xf_numpy.py b/hexrd/transforms/xf_numpy.py index 9cdb549e..d7effd59 100644 --- a/hexrd/transforms/xf_numpy.py +++ b/hexrd/transforms/xf_numpy.py @@ -27,13 +27,15 @@ # ============================================================================= # ??? do we want to set np.seterr(invalid='ignore') to avoid nan warnings? +from __future__ import absolute_import + import numpy as np from numpy import float_ as npfloat from numpy import int_ as npint # from hexrd import constants as cnst import hexrd.constants as cnst - +from .transforms_definitions import xf_api # ============================================================================= # HELPER FUNCTIONS @@ -206,17 +208,12 @@ def _z_project(x, y): # ============================================================================= +@xf_api def angles_to_gvec( angs, beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, chi=None, rmat_c=None): - """ - Takes triplets of angles in the beam frame (2*theta, eta, omega) - to components of unit G-vectors in the LAB frame. If the omega - values are not trivial (i.e. angs[:, 2] = 0.), then the components - are in the SAMPLE frame. If the crystal rmat is specified and - is not the identity, then the components are in the CRYSTAL frame. - """ + angs = np.atleast_2d(angs) nvecs, dim = angs.shape @@ -232,18 +229,12 @@ def angles_to_gvec( return _beam_to_crystal(gvec_b, beam_vec=beam_vec, eta_vec=eta_vec, rmat_s=rmat_s, rmat_c=rmat_c) - +@xf_api def angles_to_dvec( angs, beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, chi=None, rmat_c=None): - """ - Takes triplets of angles in the beam frame (2*theta, eta, omega) - to components of unit diffraction vectors in the LAB frame. If the - omega values are not trivial (i.e. angs[:, 2] = 0.), then the - components are in the SAMPLE frame. If the crystal rmat is specified - and is not the identity, then the components are in the CRYSTAL frame. - """ + angs = np.atleast_2d(angs) nvecs, dim = angs.shape @@ -260,73 +251,14 @@ def angles_to_dvec( rmat_s=rmat_s, rmat_c=rmat_c) +@xf_api def gvec_to_xy(gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, beam_vec=cnst.beam_vec, vmat_inv=None, bmat=None): - """ - Takes a concatenated list of reciprocal lattice vectors components in the - CRYSTAL FRAME to the specified detector-relative frame, subject to the - following: - - 1) it must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Parameters - ---------- - gvec_c : array_like - Concatenated triplets of G-vector components in the CRYSTAL FRAME. - rmat_d : array_like - The (3, 3) COB matrix taking components in the - DETECTOR FRAME to the LAB FRAME - rmat_s : array_like - The (3, 3) COB matrix taking components in the - SAMPLE FRAME to the LAB FRAME - rmat_c : array_like - The (3, 3) COB matrix taking components in the - CRYSTAL FRAME to the SAMPLE FRAME - tvec_d : array_like - The (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME - tvec_s : array_like - The (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME - tvec_c : array_like - The (3, ) translation vector connecting SAMPLE FRAME to CRYSTAL FRAME - beam_vec : array_like, optional - The (3, ) incident beam propagation vector components in the LAB FRAME; - the default is [0, 0, -1], which is the standard setting. - vmat_inv : array_like, optional - The (3, 3) matrix of inverse stretch tensor components in the - SAMPLE FRAME. The default is None, which implies a strain-free state - (i.e. V = I). - bmat : array_like, optional - The (3, 3) COB matrix taking components in the - RECIPROCAL LATTICE FRAME to the CRYSTAL FRAME; if supplied, it is - assumed that the input `gvecs` are G-vector components in the - RECIPROCL LATTICE FRAME (the default is None, which implies components - in the CRYSTAL FRAME) - - Returns - ------- - array_like - The (n, 2) array of [x, y] diffracted beam intersections for each of - the n input G-vectors in the DETECTOR FRAME (all Z_d coordinates are 0 - and excluded). For each input G-vector that cannot satisfy a Bragg - condition or intersect the detector plane, [NaN, Nan] is returned. - - Raises - ------ - AttributeError - The ``Raises`` section is a list of all exceptions - that are relevant to the interface. - ValueError - If `param2` is equal to `param1`. - - Notes - ----- - """ ztol = cnst.epsf # catch 1-d input case and initialize return array with NaNs @@ -396,60 +328,13 @@ def gvec_to_xy(gvec_c, return retval[:, :2] +@xf_api def xy_to_gvec(xy_d, rmat_d, rmat_s, tvec_d, tvec_s, tvec_c, rmat_b=None, distortion=None, output_ref=False): - """ - Takes a list cartesian (x, y) pairs in the DETECTOR FRAME and calculates - the associated reciprocal lattice (G) vectors and (bragg angle, azimuth) - pairs with respect to the specified beam and azimth (eta) reference - directions. - - Parameters - ---------- - xy_d : array_like - (n, 2) array of n (x, y) coordinates in DETECTOR FRAME - rmat_d : array_like - (3, 3) COB matrix taking components in the - DETECTOR FRAME to the LAB FRAME - rmat_s : array_like - (3, 3) COB matrix taking components in the - SAMPLE FRAME to the LAB FRAME - tvec_d : array_like - (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME - tvec_s : array_like - (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME - tvec_c : array_like - (3, ) translation vector connecting SAMPLE FRAME to CRYSTAL FRAME - rmat_b : array_like, optional - (3, 3) COB matrix taking components in the BEAM FRAME to the LAB FRAME; - defaults to None, which implies the standard setting of identity. - distortion : distortion class, optional - Default is None - output_ref : bool, optional - If True, prepends the apparent bragg angle and azimuth with respect to - the SAMPLE FRAME (ignoring effect of non-zero tvec_c) - - Returns - ------- - array_like - (n, 2) ndarray containing the (tth, eta) pairs associated with each - (x, y) associated with gVecs - array_like - (n, 3) ndarray containing the associated G vector directions in the - LAB FRAME - array_like, optional - if output_ref is True - - Notes - ----- - ???: is there a need to flatten the tvec inputs? - ???: include optional wavelength input for returning G with magnitude? - ???: is there a need to check that rmat_b is orthogonal if spec'd? - """ # catch 1-d input and grab number of input vectors xy_d = np.atleast_2d(xy_d) @@ -493,92 +378,9 @@ def xy_to_gvec(xy_d, return (tth, eta), ghat_l +@xf_api def solve_omega(gvecs, chi, rmat_c, wavelength, bmat=None, vmat_inv=None, rmat_b=None): - """ - For the monochromatic rotation method. - - Solve the for the rotation angle pairs that satisfy the bragg conditions - for an input list of G-vector components. - - Parameters - ---------- - gvecs : array_like - Concatenated triplets of G-vector components in either the - CRYSTAL FRAME or RECIPROCAL FRAME (see optional kwarg `bmat` below). - The shape when cast as a 2-d ndarray is (n, 3), representing n vectors. - chi : float - The inclination angle of the goniometer axis (standard coords) - rmat_c : array_like - (3, 3) COB matrix taking components in the - CRYSTAL FRAME to the SAMPLE FRAME - wavelength : float - The X-ray wavelength in Angstroms - bmat : array_like, optional - The (3, 3) COB matrix taking components in the - RECIPROCAL LATTICE FRAME to the CRYSTAL FRAME; if supplied, it is - assumed that the input `gvecs` are G-vector components in the - RECIPROCL LATTICE FRAME (the default is None, which implies components - in the CRYSTAL FRAME) - vmat_inv : array_like, optional - The (3, 3) matrix of inverse stretch tensor components in the - SAMPLE FRAME. The default is None, which implies a strain-free state - (i.e. V = I). - rmat_b : array_like, optional - (3, 3) COB matrix taking components in the BEAM FRAME to the LAB FRAME; - defaults to None, which implies the standard setting of identity. - - Returns - ------- - ome0 : array_like - The (n, 3) ndarray containing the feasible (tth, eta, ome) triplets for - each input hkl (first solution) - ome1 : array_like - The (n, 3) ndarray containing the feasible (tth, eta, ome) triplets for - each input hkl (second solution) - - Notes - ----- - The reciprocal lattice vector, G, will satisfy the the Bragg condition - when: - - b.T * G / ||G|| = -sin(theta) - - where b is the incident beam direction (k_i) and theta is the Bragg - angle consistent with G and the specified wavelength. The components of - G in the lab frame in this case are obtained using the crystal - orientation, Rc, and the single-parameter oscillation matrix, Rs(ome): - - Rs(ome) * Rc * G / ||G|| - - The equation above can be rearranged to yeild an expression of the form: - - a*sin(ome) + b*cos(ome) = c - - which is solved using the relation: - - a*sin(x) + b*cos(x) = sqrt(a**2 + b**2) * sin(x + alpha) - - --> sin(x + alpha) = c / sqrt(a**2 + b**2) - - where: - - alpha = arctan2(b, a) - - The solutions are: - - / - | arcsin(c / sqrt(a**2 + b**2)) - alpha - x = < - | pi - arcsin(c / sqrt(a**2 + b**2)) - alpha - \ - - There is a double root in the case the reflection is tangent to the - Debye-Scherrer cone (c**2 = a**2 + b**2), and no solution if the - Laue condition cannot be satisfied (filled with NaNs in the results - array here) - """ - gvecs = np.atleast_2d(gvecs) # sin and cos of the oscillation axis tilt @@ -681,13 +483,8 @@ def solve_omega(gvecs, chi, rmat_c, wavelength, # UTILITY FUNCTIONS # ============================================================================= - +@xf_api def angular_difference(ang_list0, ang_list1, units=cnst.angular_units): - """ - Do the proper (acute) angular difference in the context of a branch cut. - - *) Default angular range is [-pi, pi] - """ period = cnst.period_dict[units] # take difference as arrays diffAngles = np.atleast_1d(ang_list0) - np.atleast_1d(ang_list1) @@ -695,14 +492,8 @@ def angular_difference(ang_list0, ang_list1, units=cnst.angular_units): return abs(np.remainder(diffAngles + 0.5*period, period) - 0.5*period) +@xf_api def map_angle(ang, *args, **kwargs): - """ - Utility routine to map an angle into a specified period - - actual function is map_angle(ang[, range], units=cnst.angular_units). - range is optional and defaults to the appropriate angle for the unit - centered on 0. - """ units = cnst.angular_units period = cnst.period_dict[units] @@ -752,11 +543,8 @@ def map_angle(ang, *args, **kwargs): retval = np.mod(ang + 0.5*period, period) - 0.5*period return retval - +@xf_api def row_norm(a): - """ - normalize array of row vectors (vstacked, axis = 1) - """ if len(a.shape) > 2: raise RuntimeError( "incorrect shape: arg must be 1-d or 2-d, yours is %d" @@ -765,10 +553,8 @@ def row_norm(a): return np.sqrt(sum(np.asarray(a)**2, 1)) +@xf_api def unit_vector(a): - """ - normalize an array of row vectors (vstacked, axis=0) - """ a = np.atleast_2d(a) n = a.shape[1] @@ -779,15 +565,8 @@ def unit_vector(a): return (a/nrm).squeeze() +@xf_api def make_sample_rmat(chi, ome): - """ - Make SAMPLE frame rotation matrices as composition of - rotation of ome about the axis - - [0., cos(chi), sin(chi)] - - in the lab frame - """ # angle chi about LAB X cchi = np.cos(chi) schi = np.sin(chi) @@ -807,10 +586,8 @@ def make_sample_rmat(chi, ome): return rmat_s +@xf_api def make_rmat_of_expmap(exp_map): - """ - Calculates the rotation matrix from an exponential map - """ phi = np.sqrt( exp_map[0]*exp_map[0] + exp_map[1]*exp_map[1] @@ -831,20 +608,14 @@ def make_rmat_of_expmap(exp_map): return rmat +@xf_api def make_binary_rmat(n): - """ - make a binary rotation matrix about the specified axis - """ assert len(n) == 3, 'Axis input does not have 3 components' return 2*np.outer(n, n) - cnst.identity_3x3 +@xf_api def make_beam_rmat(bvec_l, evec_l): - """ - make eta basis COB matrix with beam antiparallel with Z - - takes components from BEAM frame to LAB - """ # normalize input bhat_l = unit_vector(bvec_l) ehat_l = unit_vector(evec_l) @@ -861,16 +632,8 @@ def make_beam_rmat(bvec_l, evec_l): return np.vstack([Xe, Ye, -bhat_l]) +@xf_api def angles_in_range(angles, starts, stops, degrees=True): - """Determine whether angles lie in or out of specified ranges - - *angles* - a list/array of angles - *starts* - a list of range starts - *stops* - a list of range stops - - OPTIONAL ARGS: - *degrees* - [True] angles & ranges in degrees (or radians) - """ tau = 360.0 if degrees else 2*np.pi nw = len(starts) na = len(angles) @@ -886,14 +649,8 @@ def angles_in_range(angles, starts, stops, degrees=True): return in_range +@xf_api def validate_angle_ranges(ang_list, startAngs, stopAngs, ccw=True): - """ - A better way to go. find out if an angle is in the range - CCW or CW from start to stop - - There is, of course, an ambigutiy if the start and stop angle are - the same; we treat them as implying 2*pi having been mapped - """ # Prefer ravel over flatten because flatten never skips the copy ang_list = np.asarray(ang_list).ravel() startAngs = np.asarray(startAngs).ravel() @@ -986,26 +743,8 @@ def validate_angle_ranges(ang_list, startAngs, stopAngs, ccw=True): return reflInRange +@xf_api def rotate_vecs_about_axis(angle, axis, vecs): - """ - Rotate vectors about an axis - - INPUTS - *angle* - array of angles (len == 1 or n) - *axis* - array of unit vectors (shape == (3, 1) or (3, n)) - *vec* - array of vectors to be rotated (shape = (3, n)) - - Quaternion formula: - if we split v into parallel and perpedicular components w.r.t. the - axis of quaternion q, - - v = a + n - - then the action of rotating the vector dot(R(q), v) becomes - - v_rot = (q0**2 - |q|**2)(a + n) + 2*dot(q, a)*q + 2*q0*cross(q, n) - - """ angle = np.atleast_1d(angle) # nvecs = vecs.shape[1] # assume column vecs @@ -1048,33 +787,8 @@ def rotate_vecs_about_axis(angle, axis, vecs): return v_rot +@xf_api def quat_product_matrix(q, mult='right'): - """ - Form 4 x 4 array to perform the quaternion product - - USAGE - qmat = quatProductMatrix(q, mult='right') - - INPUTS - 1) quats is (4,), an iterable representing a unit quaternion - horizontally concatenated - 2) mult is a keyword arg, either 'left' or 'right', denoting - the sense of the multiplication: - - / quatProductMatrix(h, mult='right') * q - q * h --> < - \ quatProductMatrix(q, mult='left') * h - - OUTPUTS - 1) qmat is (4, 4), the left or right quaternion product - operator - - NOTES - *) This function is intended to replace a cross-product based - routine for products of quaternions with large arrays of - quaternions (e.g. applying symmetries to a large set of - orientations). - """ if mult == 'right': qmat = np.array([[ q[0], -q[1], -q[2], -q[3]], [ q[1], q[0], q[3], -q[2]], @@ -1090,10 +804,8 @@ def quat_product_matrix(q, mult='right'): return qmat +@xf_api def quat_distance(q1, q2, qsym): - """ - find the distance between two unit quaternions under symmetry group - """ # qsym from PlaneData objects are (4, nsym) # convert symmetries to (4, 4) qprod matrices nsym = qsym.shape[1] From 88a837456a0e392494d30ac3fa0bc83fa75b1b6d Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 4 Jun 2018 17:33:06 +0200 Subject: [PATCH 09/30] Added skeleton test files. --- .gitignore | 1 + hexrd/transforms/tests/__init__.py | 1 + hexrd/transforms/tests/test_angle_in_range.py | 26 +++ hexrd/transforms/tests/test_angles_to_dvec.py | 26 +++ hexrd/transforms/tests/test_angles_to_gvec.py | 26 +++ .../tests/test_angular_difference.py | 26 +++ hexrd/transforms/tests/test_gvec_to_xy.py | 26 +++ hexrd/transforms/tests/test_make_beam_rmat.py | 26 +++ .../transforms/tests/test_make_binary_rmat.py | 26 +++ .../tests/test_make_rmat_of_expmap.py | 26 +++ .../transforms/tests/test_make_sample_rmat.py | 26 +++ hexrd/transforms/tests/test_map_angle.py | 26 +++ hexrd/transforms/tests/test_quat_distance.py | 26 +++ .../tests/test_quat_product_matrix.py | 26 +++ .../tests/test_rotate_vecs_about_axis.py | 26 +++ hexrd/transforms/tests/test_row_norm.py | 26 +++ hexrd/transforms/tests/test_solve_omega.py | 26 +++ hexrd/transforms/tests/test_unit_vector.py | 26 +++ .../tests/test_validate_angle_ranges.py | 26 +++ hexrd/transforms/tests/test_xy_to_gvec.py | 26 +++ hexrd/transforms/transforms_definitions.py | 11 +- hexrd/transforms/xf_capi.py | 193 ++++++++++++++++++ hexrd/transforms/xf_numba.py | 4 +- hexrd/transforms/xf_numpy.py | 2 +- 24 files changed, 675 insertions(+), 5 deletions(-) create mode 100644 hexrd/transforms/tests/__init__.py create mode 100644 hexrd/transforms/tests/test_angle_in_range.py create mode 100644 hexrd/transforms/tests/test_angles_to_dvec.py create mode 100644 hexrd/transforms/tests/test_angles_to_gvec.py create mode 100644 hexrd/transforms/tests/test_angular_difference.py create mode 100644 hexrd/transforms/tests/test_gvec_to_xy.py create mode 100644 hexrd/transforms/tests/test_make_beam_rmat.py create mode 100644 hexrd/transforms/tests/test_make_binary_rmat.py create mode 100644 hexrd/transforms/tests/test_make_rmat_of_expmap.py create mode 100644 hexrd/transforms/tests/test_make_sample_rmat.py create mode 100644 hexrd/transforms/tests/test_map_angle.py create mode 100644 hexrd/transforms/tests/test_quat_distance.py create mode 100644 hexrd/transforms/tests/test_quat_product_matrix.py create mode 100644 hexrd/transforms/tests/test_rotate_vecs_about_axis.py create mode 100644 hexrd/transforms/tests/test_row_norm.py create mode 100644 hexrd/transforms/tests/test_solve_omega.py create mode 100644 hexrd/transforms/tests/test_unit_vector.py create mode 100644 hexrd/transforms/tests/test_validate_angle_ranges.py create mode 100644 hexrd/transforms/tests/test_xy_to_gvec.py diff --git a/.gitignore b/.gitignore index 58309e9f..1e80c6f1 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ hexrd.egg-info *.*~ *# .#* +.pytest_cache \ No newline at end of file diff --git a/hexrd/transforms/tests/__init__.py b/hexrd/transforms/tests/__init__.py new file mode 100644 index 00000000..af717b81 --- /dev/null +++ b/hexrd/transforms/tests/__init__.py @@ -0,0 +1 @@ +"""Test suite for hexrd.transforms""" diff --git a/hexrd/transforms/tests/test_angle_in_range.py b/hexrd/transforms/tests/test_angle_in_range.py new file mode 100644 index 00000000..587ced03 --- /dev/null +++ b/hexrd/transforms/tests/test_angle_in_range.py @@ -0,0 +1,26 @@ +# tests for angles_in_range + +from __future__ import absolute_import + +from .. import angles_in_range as default_angles_in_range +from ..xf_numpy import angles_in_range as numpy_angles_in_range +#from ..xf_capi import angles_in_range as capi_angles_in_range +#from ..xf_numba import angles_in_range as numba_angles_in_range + +import pytest + +all_impls = pytest.mark.parametrize('angles_in_range_impl, module_name', + [(numpy_angles_in_range, 'numpy'), + #(capi_angles_in_range, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_angles_in_range, 'default')] + ) + + +@all_impls +def test_sample1(angles_in_range_impl, module_name): + pass + +@all_impls +def test_sample2(angles_in_range_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_angles_to_dvec.py b/hexrd/transforms/tests/test_angles_to_dvec.py new file mode 100644 index 00000000..a167caf4 --- /dev/null +++ b/hexrd/transforms/tests/test_angles_to_dvec.py @@ -0,0 +1,26 @@ +# tests for angles_to_dvec + +from __future__ import absolute_import + +from .. import angles_to_dvec as default_angles_to_dvec +from ..xf_numpy import angles_to_dvec as numpy_angles_to_dvec +#from ..xf_capi import angles_to_dvec as capi_angles_to_dvec +#from ..xf_numba import angles_to_dvec as numba_angles_to_dvec + +import pytest + +all_impls = pytest.mark.parametrize('angles_to_dvec_impl, module_name', + [(numpy_angles_to_dvec, 'numpy'), + #(capi_angles_to_dvec, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_angles_to_dvec, 'default')] + ) + + +@all_impls +def test_sample1(angles_to_dvec_impl, module_name): + pass + +@all_impls +def test_sample2(angles_to_dvec_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_angles_to_gvec.py b/hexrd/transforms/tests/test_angles_to_gvec.py new file mode 100644 index 00000000..ba72cd41 --- /dev/null +++ b/hexrd/transforms/tests/test_angles_to_gvec.py @@ -0,0 +1,26 @@ +# tests for angles_to_gvec + +from __future__ import absolute_import + +from .. import angles_to_gvec as default_angles_to_gvec +from ..xf_numpy import angles_to_gvec as numpy_angles_to_gvec +from ..xf_capi import angles_to_gvec as capi_angles_to_gvec +#from ..xf_numba import angles_to_gvec as numba_angles_to_gvec + +import pytest + +all_impls = pytest.mark.parametrize('angles_to_gvec_impl, module_name', + [(numpy_angles_to_gvec, 'numpy'), + (capi_angles_to_gvec, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_angles_to_gvec, 'default')] + ) + + +@all_impls +def test_sample1(angles_to_gvec_impl, module_name): + pass + +@all_impls +def test_sample2(angles_to_gvec_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_angular_difference.py b/hexrd/transforms/tests/test_angular_difference.py new file mode 100644 index 00000000..c02a1f07 --- /dev/null +++ b/hexrd/transforms/tests/test_angular_difference.py @@ -0,0 +1,26 @@ +# tests for angular_difference + +from __future__ import absolute_import + +from .. import angular_difference as default_angular_difference +from ..xf_numpy import angular_difference as numpy_angular_difference +#from ..xf_capi import angular_difference as capi_angular_difference +#from ..xf_numba import angular_difference as numba_angular_difference + +import pytest + +all_impls = pytest.mark.parametrize('angular_difference_impl, module_name', + [(numpy_angular_difference, 'numpy'), + #(capi_angular_difference, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_angular_difference, 'default')] + ) + + +@all_impls +def test_sample1(angular_difference_impl, module_name): + pass + +@all_impls +def test_sample2(angular_difference_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_gvec_to_xy.py b/hexrd/transforms/tests/test_gvec_to_xy.py new file mode 100644 index 00000000..ae5ba8ab --- /dev/null +++ b/hexrd/transforms/tests/test_gvec_to_xy.py @@ -0,0 +1,26 @@ +# tests for angles_to_gvec + +from __future__ import absolute_import + +from .. import gvec_to_xy as default_gvec_to_xy +from ..xf_numpy import gvec_to_xy as numpy_gvec_to_xy +from ..xf_capi import gvec_to_xy as capi_gvec_to_xy +#from ..xf_numba import gvec_to_xy as numba_gvec_to_xy + +import pytest + +all_impls = pytest.mark.parametrize('gvec_to_xy_impl, module_name', + [(numpy_gvec_to_xy, 'numpy'), + (capi_gvec_to_xy, 'capi'), + #(numba_gvec_to_xy, 'numba'), + (default_gvec_to_xy, 'default')] + ) + + +@all_impls +def test_sample1(gvec_to_xy_impl, module_name): + pass + +@all_impls +def test_sample2(gvec_to_xy_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_make_beam_rmat.py b/hexrd/transforms/tests/test_make_beam_rmat.py new file mode 100644 index 00000000..88ea185c --- /dev/null +++ b/hexrd/transforms/tests/test_make_beam_rmat.py @@ -0,0 +1,26 @@ +# tests for make_beam_rmat + +from __future__ import absolute_import + +from .. import make_beam_rmat as default_make_beam_rmat +from ..xf_numpy import make_beam_rmat as numpy_make_beam_rmat +#from ..xf_capi import make_beam_rmat as capi_make_beam_rmat +#from ..xf_numba import make_beam_rmat as numba_make_beam_rmat + +import pytest + +all_impls = pytest.mark.parametrize('make_beam_rmat_impl, module_name', + [(numpy_make_beam_rmat, 'numpy'), + #(capi_make_beam_rmat, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_make_beam_rmat, 'default')] + ) + + +@all_impls +def test_sample1(make_beam_rmat_impl, module_name): + pass + +@all_impls +def test_sample2(make_beam_rmat_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_make_binary_rmat.py b/hexrd/transforms/tests/test_make_binary_rmat.py new file mode 100644 index 00000000..9d685300 --- /dev/null +++ b/hexrd/transforms/tests/test_make_binary_rmat.py @@ -0,0 +1,26 @@ +# tests for make_binary_rmat + +from __future__ import absolute_import + +from .. import make_binary_rmat as default_make_binary_rmat +from ..xf_numpy import make_binary_rmat as numpy_make_binary_rmat +#from ..xf_capi import make_binary_rmat as capi_make_binary_rmat +#from ..xf_numba import make_binary_rmat as numba_make_binary_rmat + +import pytest + +all_impls = pytest.mark.parametrize('make_binary_rmat_impl, module_name', + [(numpy_make_binary_rmat, 'numpy'), + #(capi_make_binary_rmat, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_make_binary_rmat, 'default')] + ) + + +@all_impls +def test_sample1(make_binary_rmat_impl, module_name): + pass + +@all_impls +def test_sample2(make_binary_rmat_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_make_rmat_of_expmap.py b/hexrd/transforms/tests/test_make_rmat_of_expmap.py new file mode 100644 index 00000000..3dd51a58 --- /dev/null +++ b/hexrd/transforms/tests/test_make_rmat_of_expmap.py @@ -0,0 +1,26 @@ +# tests for make_rmat_of_expmap + +from __future__ import absolute_import + +from .. import make_rmat_of_expmap as default_make_rmat_of_expmap +from ..xf_numpy import make_rmat_of_expmap as numpy_make_rmat_of_expmap +#from ..xf_capi import make_rmat_of_expmap as capi_make_rmat_of_expmap +#from ..xf_numba import make_rmat_of_expmap as numba_make_rmat_of_expmap + +import pytest + +all_impls = pytest.mark.parametrize('make_rmat_of_expmap_impl, module_name', + [(numpy_make_rmat_of_expmap, 'numpy'), + #(capi_make_rmat_of_expmap, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_make_rmat_of_expmap, 'default')] + ) + + +@all_impls +def test_sample1(make_rmat_of_expmap_impl, module_name): + pass + +@all_impls +def test_sample2(make_rmat_of_expmap_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_make_sample_rmat.py b/hexrd/transforms/tests/test_make_sample_rmat.py new file mode 100644 index 00000000..8eec5688 --- /dev/null +++ b/hexrd/transforms/tests/test_make_sample_rmat.py @@ -0,0 +1,26 @@ +# tests for make_sample_rmat + +from __future__ import absolute_import + +from .. import make_sample_rmat as default_make_sample_rmat +from ..xf_numpy import make_sample_rmat as numpy_make_sample_rmat +#from ..xf_capi import make_sample_rmat as capi_make_sample_rmat +#from ..xf_numba import make_sample_rmat as numba_make_sample_rmat + +import pytest + +all_impls = pytest.mark.parametrize('make_sample_rmat_impl, module_name', + [(numpy_make_sample_rmat, 'numpy'), + #(capi_make_sample_rmat, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_make_sample_rmat, 'default')] + ) + + +@all_impls +def test_sample1(make_sample_rmat_impl, module_name): + pass + +@all_impls +def test_sample2(make_sample_rmat_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_map_angle.py b/hexrd/transforms/tests/test_map_angle.py new file mode 100644 index 00000000..a1fd859f --- /dev/null +++ b/hexrd/transforms/tests/test_map_angle.py @@ -0,0 +1,26 @@ +# tests for map_angle + +from __future__ import absolute_import + +from .. import map_angle as default_map_angle +from ..xf_numpy import map_angle as numpy_map_angle +#from ..xf_capi import map_angle as capi_map_angle +#from ..xf_numba import map_angle as numba_map_angle + +import pytest + +all_impls = pytest.mark.parametrize('map_angle_impl, module_name', + [(numpy_map_angle, 'numpy'), + #(capi_map_angle, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_map_angle, 'default')] + ) + + +@all_impls +def test_sample1(map_angle_impl, module_name): + pass + +@all_impls +def test_sample2(map_angle_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_quat_distance.py b/hexrd/transforms/tests/test_quat_distance.py new file mode 100644 index 00000000..318ce3d9 --- /dev/null +++ b/hexrd/transforms/tests/test_quat_distance.py @@ -0,0 +1,26 @@ +# tests for quat_distance + +from __future__ import absolute_import + +from .. import quat_distance as default_quat_distance +from ..xf_numpy import quat_distance as numpy_quat_distance +#from ..xf_capi import quat_distance as capi_quat_distance +#from ..xf_numba import quat_distance as numba_quat_distance + +import pytest + +all_impls = pytest.mark.parametrize('quat_distance_impl, module_name', + [(numpy_quat_distance, 'numpy'), + #(capi_quat_distance, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_quat_distance, 'default')] + ) + + +@all_impls +def test_sample1(quat_distance_impl, module_name): + pass + +@all_impls +def test_sample2(quat_distance_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_quat_product_matrix.py b/hexrd/transforms/tests/test_quat_product_matrix.py new file mode 100644 index 00000000..cf3ac0bb --- /dev/null +++ b/hexrd/transforms/tests/test_quat_product_matrix.py @@ -0,0 +1,26 @@ +# tests for quat_product_matrix + +from __future__ import absolute_import + +from .. import quat_product_matrix as default_quat_product_matrix +from ..xf_numpy import quat_product_matrix as numpy_quat_product_matrix +#from ..xf_capi import quat_product_matrix as capi_quat_product_matrix +#from ..xf_numba import quat_product_matrix as numba_quat_product_matrix + +import pytest + +all_impls = pytest.mark.parametrize('quat_product_matrix_impl, module_name', + [(numpy_quat_product_matrix, 'numpy'), + #(capi_quat_product_matrix, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_quat_product_matrix, 'default')] + ) + + +@all_impls +def test_sample1(quat_product_matrix_impl, module_name): + pass + +@all_impls +def test_sample2(quat_product_matrix_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_rotate_vecs_about_axis.py b/hexrd/transforms/tests/test_rotate_vecs_about_axis.py new file mode 100644 index 00000000..f6539ddb --- /dev/null +++ b/hexrd/transforms/tests/test_rotate_vecs_about_axis.py @@ -0,0 +1,26 @@ +# tests for rotate_vecs_about_axis + +from __future__ import absolute_import + +from .. import rotate_vecs_about_axis as default_rotate_vecs_about_axis +from ..xf_numpy import rotate_vecs_about_axis as numpy_rotate_vecs_about_axis +#from ..xf_capi import rotate_vecs_about_axis as capi_rotate_vecs_about_axis +#from ..xf_numba import rotate_vecs_about_axis as numba_rotate_vecs_about_axis + +import pytest + +all_impls = pytest.mark.parametrize('rotate_vecs_about_axis_impl, module_name', + [(numpy_rotate_vecs_about_axis, 'numpy'), + #(capi_rotate_vecs_about_axis, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_rotate_vecs_about_axis, 'default')] + ) + + +@all_impls +def test_sample1(rotate_vecs_about_axis_impl, module_name): + pass + +@all_impls +def test_sample2(rotate_vecs_about_axis_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_row_norm.py b/hexrd/transforms/tests/test_row_norm.py new file mode 100644 index 00000000..55e440b5 --- /dev/null +++ b/hexrd/transforms/tests/test_row_norm.py @@ -0,0 +1,26 @@ +# tests for row_norm + +from __future__ import absolute_import + +from .. import row_norm as default_row_norm +from ..xf_numpy import row_norm as numpy_row_norm +#from ..xf_capi import row_norm as capi_row_norm +#from ..xf_numba import row_norm as numba_row_norm + +import pytest + +all_impls = pytest.mark.parametrize('row_norm_impl, module_name', + [(numpy_row_norm, 'numpy'), + #(capi_row_norm, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_row_norm, 'default')] + ) + + +@all_impls +def test_sample1(row_norm_impl, module_name): + pass + +@all_impls +def test_sample2(row_norm_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_solve_omega.py b/hexrd/transforms/tests/test_solve_omega.py new file mode 100644 index 00000000..d23c742a --- /dev/null +++ b/hexrd/transforms/tests/test_solve_omega.py @@ -0,0 +1,26 @@ +# tests for solve_omega + +from __future__ import absolute_import + +from .. import solve_omega as default_solve_omega +from ..xf_numpy import solve_omega as numpy_solve_omega +#from ..xf_capi import solve_omega as capi_solve_omega +#from ..xf_numba import solve_omega as numba_solve_omega + +import pytest + +all_impls = pytest.mark.parametrize('solve_omega_impl, module_name', + [(numpy_solve_omega, 'numpy'), + #(capi_solve_omega, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_solve_omega, 'default')] + ) + + +@all_impls +def test_sample1(solve_omega_impl, module_name): + pass + +@all_impls +def test_sample2(solve_omega_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_unit_vector.py b/hexrd/transforms/tests/test_unit_vector.py new file mode 100644 index 00000000..706e0707 --- /dev/null +++ b/hexrd/transforms/tests/test_unit_vector.py @@ -0,0 +1,26 @@ +# tests for unit_vector + +from __future__ import absolute_import + +from .. import unit_vector as default_unit_vector +from ..xf_numpy import unit_vector as numpy_unit_vector +#from ..xf_capi import unit_vector as capi_unit_vector +#from ..xf_numba import unit_vector as numba_unit_vector + +import pytest + +all_impls = pytest.mark.parametrize('unit_vector_impl, module_name', + [(numpy_unit_vector, 'numpy'), + #(capi_unit_vector, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_unit_vector, 'default')] + ) + + +@all_impls +def test_sample1(unit_vector_impl, module_name): + pass + +@all_impls +def test_sample2(unit_vector_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_validate_angle_ranges.py b/hexrd/transforms/tests/test_validate_angle_ranges.py new file mode 100644 index 00000000..af44d217 --- /dev/null +++ b/hexrd/transforms/tests/test_validate_angle_ranges.py @@ -0,0 +1,26 @@ +# tests for validate_angle_ranges + +from __future__ import absolute_import + +from .. import validate_angle_ranges as default_validate_angle_ranges +from ..xf_numpy import validate_angle_ranges as numpy_validate_angle_ranges +#from ..xf_capi import validate_angle_ranges as capi_validate_angle_ranges +#from ..xf_numba import validate_angle_ranges as numba_validate_angle_ranges + +import pytest + +all_impls = pytest.mark.parametrize('validate_angle_ranges_impl, module_name', + [(numpy_validate_angle_ranges, 'numpy'), + #(capi_validate_angle_ranges, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_validate_angle_ranges, 'default')] + ) + + +@all_impls +def test_sample1(validate_angle_ranges_impl, module_name): + pass + +@all_impls +def test_sample2(validate_angle_ranges_impl, module_name): + pass diff --git a/hexrd/transforms/tests/test_xy_to_gvec.py b/hexrd/transforms/tests/test_xy_to_gvec.py new file mode 100644 index 00000000..b8f85598 --- /dev/null +++ b/hexrd/transforms/tests/test_xy_to_gvec.py @@ -0,0 +1,26 @@ +# tests for xy_to_gvec + +from __future__ import absolute_import + +from .. import xy_to_gvec as default_xy_to_gvec +from ..xf_numpy import xy_to_gvec as numpy_xy_to_gvec +#from ..xf_capi import xy_to_gvec as capi_xy_to_gvec +#from ..xf_numba import xy_to_gvec as numba_xy_to_gvec + +import pytest + +all_impls = pytest.mark.parametrize('xy_to_gvec_impl, module_name', + [(numpy_xy_to_gvec, 'numpy'), + #(capi_xy_to_gvec, 'capi'), + #(numba_angles_to_gvec, 'numba'), + (default_xy_to_gvec, 'default')] + ) + + +@all_impls +def test_sample1(xy_to_gvec_impl, module_name): + pass + +@all_impls +def test_sample2(xy_to_gvec_impl, module_name): + pass diff --git a/hexrd/transforms/transforms_definitions.py b/hexrd/transforms/transforms_definitions.py index 988eaf6f..81a40cd4 100644 --- a/hexrd/transforms/transforms_definitions.py +++ b/hexrd/transforms/transforms_definitions.py @@ -9,6 +9,7 @@ if they are not met. This should always be optional and incur on no overhead unless enabled, only to be used for debugging and validation purposes. """ +from __future__ import absolute_import, print_function import os import functools @@ -28,10 +29,14 @@ def _kwargs(self): """The keyword arguments the function is supposed to have""" return None - def _PRECOND(*arg, **kwarg): + @classmethod + def _PRECOND(cls, *arg, **kwarg): + print("PRECOND (", cls.__class__.__name__,")") pass - def _POSTCOND(results, *args, **kwargs): + @classmethod + def _POSTCOND(cls, results, *args, **kwargs): + print("PRECOND (", cls.__class__.__name__,")") pass @@ -523,7 +528,7 @@ def xf_api(f): # At this point use a wrapper that calls pre and post conditions if checking # is enabled, otherwise leave the function "as is". if CHECK_API: - @functools.wraps(f, assigned={__doc__: fn_def.__doc__}) + @functools.wraps(f, assigned={"__doc__": fn_def.__doc__}) def wrapper(*args, **kwargs): fn_def._PRECOND(*args, **kwargs) result = f(*args, **kwargs) diff --git a/hexrd/transforms/xf_capi.py b/hexrd/transforms/xf_capi.py index 633f8661..676a06b2 100644 --- a/hexrd/transforms/xf_capi.py +++ b/hexrd/transforms/xf_capi.py @@ -1,2 +1,195 @@ # -*- coding: utf-8 -*- +from __future__ import absolute_import, print_function, division + +from hexrd import constants as cnst +from .. import _transforms_CAPI +from .transforms_definitions import xf_api + +@xf_api +def angles_to_gvec( + angs, + beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, + chi=0., rmat_c=cnst.identity_3x3): + + angs = np.ascontiguousarray( np.atleast_2d(angs) ) + beam_vec = np.ascontiguousarray( beam_vec.flatten() ) + eta_vec = np.ascontiguousarray( eta_vec.flatten() ) + rmat_c = np.eye(3) if rmat_c is None else np.ascontiguousarray( rmat_c ) + chi = 0.0 if chi is None else float(chi) + return _transforms_CAPI.anglesToGVec(angs, + bHat_l, eHat_l, + chi, rMat_c) + + +@xf_api +def angles_to_dvec( + angs, + beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, + chi=0., rmat_c=None): + # TODO: Improve capi to avoid multiplications when rmat_c is None + angs = np.ascontiguousarray( np.atleast_2d(angs) ) + beam_vec = np.ascontiguousarray( beam_vec.flatten() ) + eta_vec = np.ascontiguousarray( eta_vec.flatten() ) + rmat_c = rmat_c if rmat_c is not None else np.ascontiguousarray( np.eye(3) ) + chi = float(chi) + + return _transforms_CAPI.anglesToDVec(angs, + bHat_l, eHat_l, + chi, rMat_c) + +#@xf_api +def makeGVector(hkl, bMat): + assert hkl.shape[0] == 3, 'hkl input must be (3, n)' + return unitVector(np.dot(bMat, hkl)) + + +@xf_api +def gvec_to_xy(gvec_c, + rmat_d, rmat_s, rmat_c, + tvec_d, tvec_s, tvec_c, + beam_vec=cnst.beam_vec, + vmat_inv=None, + bmat=None): + gvec_c = np.ascontiguousarray( np.atleast_2d( gvec_c ) ) + rmat_s = np.ascontiguousarray( rmat_s) + tvec_d = np.ascontiguousarray( tvec_d.flatten() ) + tvec_s = np.ascontiguousarray( tvec_s.flatten() ) + tvec_c = np.ascontiguousarray( tvec_c.flatten() ) + beam_vec = np.ascontiguousarray( beam_vec.flatten() ) + + # depending on the number of dimensions of rmat_s use either the array version + # or the "scalar" (over rmat_s) version. Note that rmat_s is either a 3x3 matrix + # (ndim 2) or an nx3x4 array of matrices (ndim 3) + if rmat_s.ndim > 2: + return _transforms_CAPI.gvecToDetectorXYArray(gvec_c, + rmat_d, tvec_s, tvec_c, + beam_vec) + else: + return _transforms_CAPI.gvecToDetectorXY(gvec_c, + rmat_d, rmat_s, rmat_c, + tvec_d, tvec_s, tvec_c, + beam_vec) + + +@xf_api +def xy_to_gvec(xy_d, + rmat_d, rmat_s, + tvec_d, tvec_s, tvec_c, + rmat_b=None, + distortion=None, + output_ref=False): + # This was "detectorXYToGvec" previously. There is a change in the interface where + # 'beamVec' + # beamVec? -> + # etaVec? + + # TODO: Fix this with the new interface + xy_det = np.ascontiguousarray( np.atleast_2d(xy_det) ) + tVec_d = np.ascontiguousarray( tVec_d.flatten() ) + tVec_s = np.ascontiguousarray( tVec_s.flatten() ) + tVec_c = np.ascontiguousarray( tVec_c.flatten() ) + beamVec = np.ascontiguousarray( beamVec.flatten() ) + etaVec = np.ascontiguousarray( etaVec.flatten() ) + return _transforms_CAPI.detectorXYToGvec(xy_det, + rMat_d, rMat_s, + tVec_d, tVec_s, tVec_c, + beamVec, etaVec) + + + +#@xf_api +def oscillAnglesOfHKLs(hkls, chi, rMat_c, bMat, wavelength, + vInv=None, beamVec=cnst.beam_vec, etaVec=cnst.eta_vec): + # this was oscillAnglesOfHKLs + hkls = np.array(hkls, dtype=float, order='C') + if vInv is None: + vInv = np.ascontiguousarray(vInv_ref.flatten()) + else: + vInv = np.ascontiguousarray(vInv.flatten()) + beamVec = np.ascontiguousarray(beamVec.flatten()) + etaVec = np.ascontiguousarray(etaVec.flatten()) + bMat = np.ascontiguousarray(bMat) + return _transforms_CAPI.oscillAnglesOfHKLs( + hkls, chi, rMat_c, bMat, wavelength, vInv, beamVec, etaVec + ) + + +#@xf_api +def unitRowVector(vecIn): + vecIn = np.ascontiguousarray(vecIn) + if vecIn.ndim == 1: + return _transforms_CAPI.unitRowVector(vecIn) + elif vecIn.ndim == 2: + return _transforms_CAPI.unitRowVectors(vecIn) + else: + assert vecIn.ndim in [1,2], "incorrect arg shape; must be 1-d or 2-d, yours is %d-d" % (a.ndim) + + +#@xf_api +def makeDetectorRotMat(tiltAngles): + arg = np.ascontiguousarray(np.r_[tiltAngles].flatten()) + return _transforms_CAPI.makeDetectorRotMat(arg) + + +#@xf_api +def makeOscillRotMat(oscillAngles): + arg = np.ascontiguousarray(np.r_[oscillAngles].flatten()) + return _transforms_CAPI.makeOscillRotMat(arg) + + +#@xf_api +def makeOscillRotMatArray(chi, omeArray): + arg = np.ascontiguousarray(omeArray) + return _transforms_CAPI.makeOscillRotMatArray(chi, arg) + + +#@xf_api +def makeRotMatOfExpMap(expMap): + arg = np.ascontiguousarray(expMap.flatten()) + return _transforms_CAPI.makeRotMatOfExpMap(arg) + + +#@xf_api +def makeRotMatOfQuat(quats): + arg = np.ascontiguousarray(quats) + return _transforms_CAPI.makeRotMatOfQuat(arg) + + +#@xf_api +def makeBinaryRotMat(axis): + arg = np.ascontiguousarray(axis.flatten()) + return _transforms_CAPI.makeBinaryRotMat(arg) + + +#@xf_api +def makeEtaFrameRotMat(bHat_l, eHat_l): + arg1 = np.ascontiguousarray(bHat_l.flatten()) + arg2 = np.ascontiguousarray(eHat_l.flatten()) + return _transforms_CAPI.makeEtaFrameRotMat(arg1, arg2) + + +#@xf_api +def validateAngleRanges(angList, angMin, angMax, ccw=True): + angList = angList.astype(np.double, order="C") + angMin = angMin.astype(np.double, order="C") + angMax = angMax.astype(np.double, order="C") + return _transforms_CAPI.validateAngleRanges(angList,angMin,angMax,ccw) + + +#@xf_api +def rotate_vecs_about_axis(angle, axis, vecs): + return _transforms_CAPI.rotate_vecs_about_axis(angle, axis, vecs) + + +#@xf_api +def quat_distance(q1, q2, qsym): + q1 = np.ascontiguousarray(q1.flatten()) + q2 = np.ascontiguousarray(q2.flatten()) + return _transforms_CAPI.quat_distance(q1, q2, qsym) + + +#@xf_api +def homochoricOfQuat(quats): + q = np.ascontiguousarray(quats.T) + return _transforms_CAPI.homochoricOfQuat(q) diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index 38d0b386..3dbbb820 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -27,12 +27,14 @@ # ============================================================================= # ??? do we want to set np.seterr(invalid='ignore') to avoid nan warnings? +from __future__ import absolute_import + import numpy as np from numpy import float_ as npfloat from numpy import int_ as npint # from hexrd import constants as cnst -import constants as cnst +from .. import constants as cnst import numba diff --git a/hexrd/transforms/xf_numpy.py b/hexrd/transforms/xf_numpy.py index d7effd59..61a90f4d 100644 --- a/hexrd/transforms/xf_numpy.py +++ b/hexrd/transforms/xf_numpy.py @@ -27,7 +27,7 @@ # ============================================================================= # ??? do we want to set np.seterr(invalid='ignore') to avoid nan warnings? -from __future__ import absolute_import +from __future__ import absolute_import, print_function, division import numpy as np from numpy import float_ as npfloat From b4bfe16fa38d561c025b625706d559c8d348887b Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 4 Jun 2018 17:48:54 +0200 Subject: [PATCH 10/30] explained a bit the transforms testing. Signed-off-by: Oscar Villellas --- hexrd/transforms/tests/__init__.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/hexrd/transforms/tests/__init__.py b/hexrd/transforms/tests/__init__.py index af717b81..5e7bd357 100644 --- a/hexrd/transforms/tests/__init__.py +++ b/hexrd/transforms/tests/__init__.py @@ -1 +1,17 @@ """Test suite for hexrd.transforms""" + +# Use a file for each api function to test named test_.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. + + + From 996a1fc9211a5944ef017648f00e4a9683700295 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 15 Oct 2018 17:52:09 +0200 Subject: [PATCH 11/30] Added signature checking to API. Performed a first pass on xf_capi. --- hexrd/transforms/tests/test_angles_to_dvec.py | 4 +- hexrd/transforms/tests/test_make_beam_rmat.py | 4 +- .../transforms/tests/test_make_binary_rmat.py | 4 +- .../tests/test_make_rmat_of_expmap.py | 4 +- hexrd/transforms/tests/test_quat_distance.py | 4 +- .../tests/test_rotate_vecs_about_axis.py | 4 +- hexrd/transforms/tests/test_unit_vector.py | 4 +- .../tests/test_validate_angle_ranges.py | 4 +- hexrd/transforms/tests/test_xy_to_gvec.py | 4 +- hexrd/transforms/transforms_definitions.py | 135 +++++++++++++++- hexrd/transforms/xf_capi.py | 153 ++++++++++++------ hexrd/transforms/xf_numpy.py | 68 ++++---- 12 files changed, 281 insertions(+), 111 deletions(-) diff --git a/hexrd/transforms/tests/test_angles_to_dvec.py b/hexrd/transforms/tests/test_angles_to_dvec.py index a167caf4..f4d32432 100644 --- a/hexrd/transforms/tests/test_angles_to_dvec.py +++ b/hexrd/transforms/tests/test_angles_to_dvec.py @@ -4,14 +4,14 @@ from .. import angles_to_dvec as default_angles_to_dvec from ..xf_numpy import angles_to_dvec as numpy_angles_to_dvec -#from ..xf_capi import angles_to_dvec as capi_angles_to_dvec +from ..xf_capi import angles_to_dvec as capi_angles_to_dvec #from ..xf_numba import angles_to_dvec as numba_angles_to_dvec import pytest all_impls = pytest.mark.parametrize('angles_to_dvec_impl, module_name', [(numpy_angles_to_dvec, 'numpy'), - #(capi_angles_to_dvec, 'capi'), + (capi_angles_to_dvec, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_angles_to_dvec, 'default')] ) diff --git a/hexrd/transforms/tests/test_make_beam_rmat.py b/hexrd/transforms/tests/test_make_beam_rmat.py index 88ea185c..95b9522a 100644 --- a/hexrd/transforms/tests/test_make_beam_rmat.py +++ b/hexrd/transforms/tests/test_make_beam_rmat.py @@ -4,14 +4,14 @@ from .. import make_beam_rmat as default_make_beam_rmat from ..xf_numpy import make_beam_rmat as numpy_make_beam_rmat -#from ..xf_capi import make_beam_rmat as capi_make_beam_rmat +from ..xf_capi import make_beam_rmat as capi_make_beam_rmat #from ..xf_numba import make_beam_rmat as numba_make_beam_rmat import pytest all_impls = pytest.mark.parametrize('make_beam_rmat_impl, module_name', [(numpy_make_beam_rmat, 'numpy'), - #(capi_make_beam_rmat, 'capi'), + (capi_make_beam_rmat, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_make_beam_rmat, 'default')] ) diff --git a/hexrd/transforms/tests/test_make_binary_rmat.py b/hexrd/transforms/tests/test_make_binary_rmat.py index 9d685300..4b5be96b 100644 --- a/hexrd/transforms/tests/test_make_binary_rmat.py +++ b/hexrd/transforms/tests/test_make_binary_rmat.py @@ -4,14 +4,14 @@ from .. import make_binary_rmat as default_make_binary_rmat from ..xf_numpy import make_binary_rmat as numpy_make_binary_rmat -#from ..xf_capi import make_binary_rmat as capi_make_binary_rmat +from ..xf_capi import make_binary_rmat as capi_make_binary_rmat #from ..xf_numba import make_binary_rmat as numba_make_binary_rmat import pytest all_impls = pytest.mark.parametrize('make_binary_rmat_impl, module_name', [(numpy_make_binary_rmat, 'numpy'), - #(capi_make_binary_rmat, 'capi'), + (capi_make_binary_rmat, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_make_binary_rmat, 'default')] ) diff --git a/hexrd/transforms/tests/test_make_rmat_of_expmap.py b/hexrd/transforms/tests/test_make_rmat_of_expmap.py index 3dd51a58..ac6b4074 100644 --- a/hexrd/transforms/tests/test_make_rmat_of_expmap.py +++ b/hexrd/transforms/tests/test_make_rmat_of_expmap.py @@ -4,14 +4,14 @@ from .. import make_rmat_of_expmap as default_make_rmat_of_expmap from ..xf_numpy import make_rmat_of_expmap as numpy_make_rmat_of_expmap -#from ..xf_capi import make_rmat_of_expmap as capi_make_rmat_of_expmap +from ..xf_capi import make_rmat_of_expmap as capi_make_rmat_of_expmap #from ..xf_numba import make_rmat_of_expmap as numba_make_rmat_of_expmap import pytest all_impls = pytest.mark.parametrize('make_rmat_of_expmap_impl, module_name', [(numpy_make_rmat_of_expmap, 'numpy'), - #(capi_make_rmat_of_expmap, 'capi'), + (capi_make_rmat_of_expmap, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_make_rmat_of_expmap, 'default')] ) diff --git a/hexrd/transforms/tests/test_quat_distance.py b/hexrd/transforms/tests/test_quat_distance.py index 318ce3d9..3a8ec97c 100644 --- a/hexrd/transforms/tests/test_quat_distance.py +++ b/hexrd/transforms/tests/test_quat_distance.py @@ -4,14 +4,14 @@ from .. import quat_distance as default_quat_distance from ..xf_numpy import quat_distance as numpy_quat_distance -#from ..xf_capi import quat_distance as capi_quat_distance +from ..xf_capi import quat_distance as capi_quat_distance #from ..xf_numba import quat_distance as numba_quat_distance import pytest all_impls = pytest.mark.parametrize('quat_distance_impl, module_name', [(numpy_quat_distance, 'numpy'), - #(capi_quat_distance, 'capi'), + (capi_quat_distance, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_quat_distance, 'default')] ) diff --git a/hexrd/transforms/tests/test_rotate_vecs_about_axis.py b/hexrd/transforms/tests/test_rotate_vecs_about_axis.py index f6539ddb..af341c93 100644 --- a/hexrd/transforms/tests/test_rotate_vecs_about_axis.py +++ b/hexrd/transforms/tests/test_rotate_vecs_about_axis.py @@ -4,14 +4,14 @@ from .. import rotate_vecs_about_axis as default_rotate_vecs_about_axis from ..xf_numpy import rotate_vecs_about_axis as numpy_rotate_vecs_about_axis -#from ..xf_capi import rotate_vecs_about_axis as capi_rotate_vecs_about_axis +from ..xf_capi import rotate_vecs_about_axis as capi_rotate_vecs_about_axis #from ..xf_numba import rotate_vecs_about_axis as numba_rotate_vecs_about_axis import pytest all_impls = pytest.mark.parametrize('rotate_vecs_about_axis_impl, module_name', [(numpy_rotate_vecs_about_axis, 'numpy'), - #(capi_rotate_vecs_about_axis, 'capi'), + (capi_rotate_vecs_about_axis, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_rotate_vecs_about_axis, 'default')] ) diff --git a/hexrd/transforms/tests/test_unit_vector.py b/hexrd/transforms/tests/test_unit_vector.py index 706e0707..bf140840 100644 --- a/hexrd/transforms/tests/test_unit_vector.py +++ b/hexrd/transforms/tests/test_unit_vector.py @@ -4,14 +4,14 @@ from .. import unit_vector as default_unit_vector from ..xf_numpy import unit_vector as numpy_unit_vector -#from ..xf_capi import unit_vector as capi_unit_vector +from ..xf_capi import unit_vector as capi_unit_vector #from ..xf_numba import unit_vector as numba_unit_vector import pytest all_impls = pytest.mark.parametrize('unit_vector_impl, module_name', [(numpy_unit_vector, 'numpy'), - #(capi_unit_vector, 'capi'), + (capi_unit_vector, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_unit_vector, 'default')] ) diff --git a/hexrd/transforms/tests/test_validate_angle_ranges.py b/hexrd/transforms/tests/test_validate_angle_ranges.py index af44d217..14f2084a 100644 --- a/hexrd/transforms/tests/test_validate_angle_ranges.py +++ b/hexrd/transforms/tests/test_validate_angle_ranges.py @@ -4,14 +4,14 @@ from .. import validate_angle_ranges as default_validate_angle_ranges from ..xf_numpy import validate_angle_ranges as numpy_validate_angle_ranges -#from ..xf_capi import validate_angle_ranges as capi_validate_angle_ranges +from ..xf_capi import validate_angle_ranges as capi_validate_angle_ranges #from ..xf_numba import validate_angle_ranges as numba_validate_angle_ranges import pytest all_impls = pytest.mark.parametrize('validate_angle_ranges_impl, module_name', [(numpy_validate_angle_ranges, 'numpy'), - #(capi_validate_angle_ranges, 'capi'), + (capi_validate_angle_ranges, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_validate_angle_ranges, 'default')] ) diff --git a/hexrd/transforms/tests/test_xy_to_gvec.py b/hexrd/transforms/tests/test_xy_to_gvec.py index b8f85598..a2fdb89a 100644 --- a/hexrd/transforms/tests/test_xy_to_gvec.py +++ b/hexrd/transforms/tests/test_xy_to_gvec.py @@ -4,14 +4,14 @@ from .. import xy_to_gvec as default_xy_to_gvec from ..xf_numpy import xy_to_gvec as numpy_xy_to_gvec -#from ..xf_capi import xy_to_gvec as capi_xy_to_gvec +from ..xf_capi import xy_to_gvec as capi_xy_to_gvec #from ..xf_numba import xy_to_gvec as numba_xy_to_gvec import pytest all_impls = pytest.mark.parametrize('xy_to_gvec_impl, module_name', [(numpy_xy_to_gvec, 'numpy'), - #(capi_xy_to_gvec, 'capi'), + (capi_xy_to_gvec, 'capi'), #(numba_angles_to_gvec, 'numba'), (default_xy_to_gvec, 'default')] ) diff --git a/hexrd/transforms/transforms_definitions.py b/hexrd/transforms/transforms_definitions.py index 81a40cd4..3535d9b3 100644 --- a/hexrd/transforms/transforms_definitions.py +++ b/hexrd/transforms/transforms_definitions.py @@ -8,17 +8,41 @@ execute the scripts forcing those conditions to be evaluated and raise errors if they are not met. This should always be optional and incur on no overhead unless enabled, only to be used for debugging and validation purposes. + +Checking of signature definitions has been added. This happens if CHECK_API is +enabled (via the HEXRD_XF_CHECK environment variable). This is implemented via +the signature method in inspect (Python 3.3+). If not available, it falls back +to the backport in funcsigs. If unavailable, CHECK_API is disabled. """ from __future__ import absolute_import, print_function import os import functools +from .. import constants as cnst + CHECK_API = os.getenv("HEXRD_XF_CHECK") +try: + from inspect import signature as get_signature +except ImportError: + try: + from funcsigs import signature as get_signature + except: + import warnings + + warnings.warn("Failed to import from inspect/funcsigs." + "Transforms API signature checking disabled.") + get_signature = None + class DEF_Func(object): """Documentation to use for the function""" + def _signature(): + """The signature of this method defines the one for the API + including default values.""" + pass + @property def _args(self): """The position arguments the function is supposed to have""" @@ -51,7 +75,17 @@ class DEF_angles_to_gvec(DEF_Func): values are not trivial (i.e. angs[:, 2] = 0.), then the components are in the SAMPLE frame. If the crystal rmat is specified and is not the identity, then the components are in the CRYSTAL frame. + + default beam_vec is defined in hexrd.constants.beam_vec + default eta_vec is defined in hexrd.constants.eta_vec """ + def _signature(angs, + beam_vec=None, + eta_vec=None, + chi=None, + rmat_c=None): + pass + @property def _args(self): return ('angs',) @@ -68,7 +102,17 @@ class DEF_angles_to_dvec(DEF_Func): omega values are not trivial (i.e. angs[:, 2] = 0.), then the components are in the SAMPLE frame. If the crystal rmat is specified and is not the identity, then the components are in the CRYSTAL frame. + + default beam_vec is defined in hexrd.constants.beam_vec + default eta_vec is defined in hexrd.constants.eta_vec """ + def _signature(angs, + beam_vec=None, + eta_vec=None, + chi=None, + rmat_c=None): + pass + @property def _args(self): return ('angs',) @@ -140,6 +184,14 @@ class DEF_gvec_to_xy(DEF_Func): ----- """ + def _signature(gvec_c, + rmat_d, rmat_s, rmat_c, + tvec_d, tvec_s, tvec_c, + beam_vec=None, + vmat_inv=None, + bmat=None): + pass + @property def _args(self): return ('gvec_c', 'rmat_d', 'rmat_s', 'rmat_c', 'tvec_d', 'tvec_s', @@ -199,6 +251,13 @@ class DEF_xy_to_gvec(DEF_Func): ???: include optional wavelength input for returning G with magnitude? ???: is there a need to check that rmat_b is orthogonal if spec'd? """ + def _signature(xy_d, + rmat_d, rmat_s, + tvec_d, tvec_s, tvec_c, + rmat_b=None, + distortion=None, + output_ref=False): + pass @property def _args(self): @@ -293,6 +352,10 @@ class DEF_solve_omega(DEF_Func): Laue condition cannot be satisfied (filled with NaNs in the results array here) """ + def _signature(gvecs, chi, rmat_c, wavelength, + bmat=None, vmat_inv=None, rmat_b=None): + pass + @property def _args(self): return ('gvecs', 'chi', 'rmat_c', 'wavelength') @@ -312,6 +375,9 @@ class DEF_angular_difference(DEF_Func): *) Default angular range is [-pi, pi] """ + def _signature(ang_list0, ang_list1, units=cnst.angular_units): + pass + @property def _args(self): return ('ang_list0', 'ang_list1') @@ -328,7 +394,13 @@ class DEF_map_angle(DEF_Func): actual function is map_angle(ang[, range], units=cnst.angular_units). range is optional and defaults to the appropriate angle for the unit centered on 0. + + accepted units are: 'radians' and 'degrees' """ + + def _signature(ang, range=None, units=cnst.angular_units): + pass + @property def _args(self): return ('ang',) @@ -342,6 +414,9 @@ class DEF_row_norm(DEF_Func): """ normalize array of row vectors (vstacked, axis = 1) """ + def _signature(a): + pass + @property def _args(self): return ('a') @@ -351,9 +426,12 @@ class DEF_unit_vector(DEF_Func): """ normalize an array of row vectors (vstacked, axis=0) """ + def _signature(vec_in): + pass + @property def _args(self): - return('a') + return('vec_in') class DEF_make_sample_rmat(DEF_Func): @@ -365,6 +443,9 @@ class DEF_make_sample_rmat(DEF_Func): in the lab frame """ + def _signature(chi, ome): + pass + @property def _args(self): return ('chi', 'ome') @@ -374,6 +455,9 @@ class DEF_make_rmat_of_expmap(DEF_Func): """ Calculates the rotation matrix from an exponential map """ + def _signature(exp_map): + pass + @property def _args(self): return ('exp_map',) @@ -383,9 +467,12 @@ class DEF_make_binary_rmat(DEF_Func): """ make a binary rotation matrix about the specified axis """ + def _signature(axis): + pass + @property def _args(self): - return ('n',) + return ('axis',) class DEF_make_beam_rmat(DEF_Func): @@ -394,6 +481,9 @@ class DEF_make_beam_rmat(DEF_Func): takes components from BEAM frame to LAB """ + def _signature(bvec_l, evec_l): + pass + @property def _args(self): return ('bvec_l', 'evec_l') @@ -409,6 +499,9 @@ class DEF_angles_in_range(DEF_Func): OPTIONAL ARGS: *degrees* - [True] angles & ranges in degrees (or radians) """ + def _signature(angles, starts, stops, degrees=True): + pass + @property def _args(self): return ('angles', 'starts', 'stops') @@ -425,9 +518,12 @@ class DEF_validate_angle_ranges(DEF_Func): There is, of course, an ambigutiy if the start and stop angle are the same; we treat them as implying 2*pi having been mapped """ + def _signature(ang_list, start_angs, stop_angs, ccw=True): + pass + @property def _args(self): - return ('ang_list', 'startAngs', 'stopAngs') + return ('ang_list', 'start_angs', 'stop_angs') @property def _kwargs(self): @@ -454,6 +550,9 @@ class DEF_rotate_vecs_about_axis(DEF_Func): v_rot = (q0**2 - |q|**2)(a + n) + 2*dot(q, a)*q + 2*q0*cross(q, n) """ + def _signature(angle, axis, vecs): + pass + @property def _args(self): return ('angle', 'axis', 'vecs') @@ -486,6 +585,9 @@ class DEF_quat_product_matrix(DEF_Func): quaternions (e.g. applying symmetries to a large set of orientations). """ + def _signature(q, mult='right'): + pass + @property def _args(self): return ('q',) @@ -499,6 +601,9 @@ class DEF_quat_distance(DEF_Func): """ find the distance between two unit quaternions under symmetry group """ + def _signature(q1, q2, qsym): + pass + @property def _args(self): return ('q1', 'q2', 'qsym') @@ -515,15 +620,33 @@ def xf_api(f): try: fn_def = globals()['DEF_'+api_call] except KeyError: - raise RuntimeError("xf_api function '%s' doesn't have a definition.") + # This happens if there is no definition for the decorated function + raise RuntimeError("'%s' definition not found." % api_call) try: if not (isinstance(fn_def.__doc__, basestring) and callable(fn_def._PRECOND) and - callable(fn_def._POSTCOND)): + callable(fn_def._POSTCOND) and + callable(fn_def._signature)): raise Exception() except Exception: - raise RuntimeError("xf_api definition for function '%s' seems incorrect.") + # A valid definition requires a string doc, and callable _PRECOND, + # _POSTCOND and _signature. + # + # __doc__ will become the decorated function's documentation. + # _PRECOND will be run on every call with args and kwargs + # _POSTCOND will be run on every call with result, args and kwargs + # _signature will be used to enforce a signature on implementations. + # + # _PRECOND and _POSTCOND will only be called if CHECK_API is enabled, + # as they will slow down execution. + raise RuntimeError("'%s' definition error." % api_call) + + # Sanity check: make sure the decorated function has the expected signature. + if get_signature is not None: + # Check that the function has the right signature + if get_signature(fn_def._signature) != get_signature(f): + raise RuntimeError("'%s' signature mismatch." % api_call) # At this point use a wrapper that calls pre and post conditions if checking # is enabled, otherwise leave the function "as is". diff --git a/hexrd/transforms/xf_capi.py b/hexrd/transforms/xf_capi.py index 676a06b2..b08e7100 100644 --- a/hexrd/transforms/xf_capi.py +++ b/hexrd/transforms/xf_capi.py @@ -1,5 +1,27 @@ # -*- coding: utf-8 -*- - +""" +Transforms module implementation using a support C extension module. + +Currently, this implementation contains code for the following functions: + + - angles_to_gvec + - angles_to_dvec + - gvec_to_xy + - xy_to_gvec (partial) + + - unit_vector + - make_rmat_of_exp_map + - make_binary_rmat + - make_beam_rmat + - validate_angle_ranges + - rotate_vecs_about_axis + - quat_distance + +There are also some functions that maybe would be needed in the transforms module: + - makeGVector + - makeRotMatOfQuat + - homochoricOfQuat +""" from __future__ import absolute_import, print_function, division from hexrd import constants as cnst @@ -9,36 +31,43 @@ @xf_api def angles_to_gvec( angs, - beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, - chi=0., rmat_c=cnst.identity_3x3): + beam_vec=None, eta_vec=None, + chi=None, rmat_c=None): + # TODO: revise + + beam_vec = beam_vec if beam_vec is not None else const.beam_vec + eta_vec = eta_vec if eta_vec is not None else const.eta_vec angs = np.ascontiguousarray( np.atleast_2d(angs) ) beam_vec = np.ascontiguousarray( beam_vec.flatten() ) eta_vec = np.ascontiguousarray( eta_vec.flatten() ) - rmat_c = np.eye(3) if rmat_c is None else np.ascontiguousarray( rmat_c ) + rmat_c = const.identity_3x3 if rmat_c is None else np.ascontiguousarray( rmat_c ) chi = 0.0 if chi is None else float(chi) + return _transforms_CAPI.anglesToGVec(angs, - bHat_l, eHat_l, - chi, rMat_c) + beam_vec, eta_vec, + chi, rmat_c) @xf_api def angles_to_dvec( angs, - beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, - chi=0., rmat_c=None): + beam_vec=None, eta_vec=None, + chi=None, rmat_c=None): # TODO: Improve capi to avoid multiplications when rmat_c is None + beam_vec = beam_vec if beam_vec is not None else const.beam_vec + eta_vec = eta_vec if eta_vec is not None else const.eta_vec + angs = np.ascontiguousarray( np.atleast_2d(angs) ) beam_vec = np.ascontiguousarray( beam_vec.flatten() ) eta_vec = np.ascontiguousarray( eta_vec.flatten() ) - rmat_c = rmat_c if rmat_c is not None else np.ascontiguousarray( np.eye(3) ) - chi = float(chi) + rmat_c = np.ascontiguousarray(rmat_c) if rmat_c is not None else const.identity_3x3 + chi = 0.0 if chi is None else float(chi) return _transforms_CAPI.anglesToDVec(angs, - bHat_l, eHat_l, - chi, rMat_c) + beam_vec, eta_vec, + chi, rmat_c) -#@xf_api def makeGVector(hkl, bMat): assert hkl.shape[0] == 3, 'hkl input must be (3, n)' return unitVector(np.dot(bMat, hkl)) @@ -48,11 +77,13 @@ def makeGVector(hkl, bMat): def gvec_to_xy(gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, - beam_vec=cnst.beam_vec, + beam_vec=None, vmat_inv=None, bmat=None): + beam_vec = beam_vec if beam_vec is not None else const.beam_vec + gvec_c = np.ascontiguousarray( np.atleast_2d( gvec_c ) ) - rmat_s = np.ascontiguousarray( rmat_s) + rmat_s = np.ascontiguousarray( rmat_s ) tvec_d = np.ascontiguousarray( tvec_d.flatten() ) tvec_s = np.ascontiguousarray( tvec_s.flatten() ) tvec_c = np.ascontiguousarray( tvec_c.flatten() ) @@ -79,22 +110,32 @@ def xy_to_gvec(xy_d, rmat_b=None, distortion=None, output_ref=False): - # This was "detectorXYToGvec" previously. There is a change in the interface where - # 'beamVec' - # beamVec? -> - # etaVec? - - # TODO: Fix this with the new interface - xy_det = np.ascontiguousarray( np.atleast_2d(xy_det) ) - tVec_d = np.ascontiguousarray( tVec_d.flatten() ) - tVec_s = np.ascontiguousarray( tVec_s.flatten() ) - tVec_c = np.ascontiguousarray( tVec_c.flatten() ) - beamVec = np.ascontiguousarray( beamVec.flatten() ) - etaVec = np.ascontiguousarray( etaVec.flatten() ) + # in the C library beam vector and eta vector are expected. However we receive + # rmat_b. Please check this! + # + # It also seems that the output_ref version is not present as the argument gets + # ignored + + rmat_b = rmat_b if rmat_b is not None else cnst.identity_3x3 + + # the code seems to ignore this argument, assume output_ref == True not implemented + assert not output_ref + + if distortion is not None: + xy_d = distortion.unwarp(xy_d) + + xy_d = np.ascontiguousarray( np.atleast_2d(xy_d) ) + rmat_d = np.ascontiguousarray( rmat_d ) + rmat_s = np.ascontiguousarray( rmat_s ) + tvec_d = np.ascontiguousarray( tvec_d.flatten() ) + tvec_s = np.ascontiguousarray( tvec_s.flatten() ) + tvec_c = np.ascontiguousarray( tvec_c.flatten() ) + beam_vec = np.ascontiguousarray( (-rmat_b[:,2]).flatten() ) + eta_vec = np.ascontiguousarray( rmat_b[:,0].flatten() ) #check this! return _transforms_CAPI.detectorXYToGvec(xy_det, - rMat_d, rMat_s, - tVec_d, tVec_s, tVec_c, - beamVec, etaVec) + rmat_d, rmat_s, + tvec_d, tvec_s, tvec_c, + beam_vec, eta_vec) @@ -115,15 +156,15 @@ def oscillAnglesOfHKLs(hkls, chi, rMat_c, bMat, wavelength, ) -#@xf_api -def unitRowVector(vecIn): - vecIn = np.ascontiguousarray(vecIn) - if vecIn.ndim == 1: - return _transforms_CAPI.unitRowVector(vecIn) - elif vecIn.ndim == 2: - return _transforms_CAPI.unitRowVectors(vecIn) +@xf_api +def unit_vector(vec_in): + vec_in = np.ascontiguousarray(vec_in) + if vec_in.ndim == 1: + return _transforms_CAPI.unitRowVector(vec_in) + elif vec_in.ndim == 2: + return _transforms_CAPI.unitRowVectors(vec_in) else: - assert vecIn.ndim in [1,2], "incorrect arg shape; must be 1-d or 2-d, yours is %d-d" % (a.ndim) + assert vec_in.ndim in [1,2], "incorrect arg shape; must be 1-d or 2-d, yours is %d-d" % (a.ndim) #@xf_api @@ -132,6 +173,9 @@ def makeDetectorRotMat(tiltAngles): return _transforms_CAPI.makeDetectorRotMat(arg) +# make_sample_rmat in CAPI is split between makeOscillRotMat +# and makeOscillRotMatArray... + #@xf_api def makeOscillRotMat(oscillAngles): arg = np.ascontiguousarray(np.r_[oscillAngles].flatten()) @@ -144,45 +188,48 @@ def makeOscillRotMatArray(chi, omeArray): return _transforms_CAPI.makeOscillRotMatArray(chi, arg) -#@xf_api -def makeRotMatOfExpMap(expMap): - arg = np.ascontiguousarray(expMap.flatten()) +@xf_api +def make_rmat_of_expmap(exp_map): + arg = np.ascontiguousarray(exp_map.flatten()) return _transforms_CAPI.makeRotMatOfExpMap(arg) +# Not in transforms API #@xf_api def makeRotMatOfQuat(quats): arg = np.ascontiguousarray(quats) return _transforms_CAPI.makeRotMatOfQuat(arg) -#@xf_api -def makeBinaryRotMat(axis): +@xf_api +def make_binary_rmat(axis): arg = np.ascontiguousarray(axis.flatten()) return _transforms_CAPI.makeBinaryRotMat(arg) -#@xf_api -def makeEtaFrameRotMat(bHat_l, eHat_l): +@xf_api +def make_beam_rmat(bvec_l, evec_l): arg1 = np.ascontiguousarray(bHat_l.flatten()) arg2 = np.ascontiguousarray(eHat_l.flatten()) return _transforms_CAPI.makeEtaFrameRotMat(arg1, arg2) -#@xf_api -def validateAngleRanges(angList, angMin, angMax, ccw=True): - angList = angList.astype(np.double, order="C") - angMin = angMin.astype(np.double, order="C") - angMax = angMax.astype(np.double, order="C") - return _transforms_CAPI.validateAngleRanges(angList,angMin,angMax,ccw) +@xf_api +def validate_angle_ranges(ang_list, start_angs, stop_angs, ccw=True): + ang_list = ang_list.astype(np.double, order="C") + start_angs = start_angs.astype(np.double, order="C") + stop_angs = stop_angs.astype(np.double, order="C") + return _transforms_CAPI.validateAngleRanges(ang_list, start_angs, stop_angs, ccw) -#@xf_api + +@xf_api def rotate_vecs_about_axis(angle, axis, vecs): + # TODO: arguments may need preparation. return _transforms_CAPI.rotate_vecs_about_axis(angle, axis, vecs) -#@xf_api +@xf_api def quat_distance(q1, q2, qsym): q1 = np.ascontiguousarray(q1.flatten()) q2 = np.ascontiguousarray(q2.flatten()) diff --git a/hexrd/transforms/xf_numpy.py b/hexrd/transforms/xf_numpy.py index 51807b4b..e1b1db79 100644 --- a/hexrd/transforms/xf_numpy.py +++ b/hexrd/transforms/xf_numpy.py @@ -33,8 +33,7 @@ from numpy import float_ as npfloat from numpy import int_ as npint -# from hexrd import constants as cnst -import hexrd.constants as cnst +from .. import constants as cnst from .transforms_definitions import xf_api # ============================================================================= @@ -211,9 +210,12 @@ def _z_project(x, y): @xf_api def angles_to_gvec( angs, - beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, + beam_vec=None, eta_vec=None, chi=None, rmat_c=None): + beam_vec = beam_vec if beam_vec is not None else const.beam_vec + eta_vec = eta_vec if eta_vec is not None else const.eta_vec + angs = np.atleast_2d(angs) nvecs, dim = angs.shape @@ -236,9 +238,12 @@ def angles_to_gvec( @xf_api def angles_to_dvec( angs, - beam_vec=cnst.beam_vec, eta_vec=cnst.eta_vec, + beam_vec=None, eta_vec=None, chi=None, rmat_c=None): + beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec + eta_vec = eta_vec if eta_vec is not None else cnst.eta_vec + angs = np.atleast_2d(angs) nvecs, dim = angs.shape @@ -263,10 +268,12 @@ def angles_to_dvec( def gvec_to_xy(gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, - beam_vec=cnst.beam_vec, + beam_vec=None, vmat_inv=None, bmat=None): + beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec + ztol = cnst.epsf # catch 1-d input case and initialize return array with NaNs @@ -503,30 +510,18 @@ def angular_difference(ang_list0, ang_list1, units=cnst.angular_units): @xf_api -def map_angle(ang, *args, **kwargs): - units = cnst.angular_units - period = cnst.period_dict[units] - - kwargKeys = kwargs.keys() - for iArg in range(len(kwargKeys)): - if kwargKeys[iArg] == 'units': - units = kwargs[kwargKeys[iArg]] - else: - raise RuntimeError( - "Unknown keyword argument: " + str(kwargKeys[iArg]) - ) +def map_angle(ang, range=None, units=cnst.angular_units): try: - period = cnst.period_dict[units.lower()] - except(KeyError): - raise RuntimeError( - "unknown angular units: " + str(kwargs[kwargKeys[iArg]]) - ) + period = cnst.period_dict(units) + except KeyError: + raise ValueError('Accepted units: %s'.format(','.join(cnst.period_dict.keys()))) + ang = np.atleast_1d(npfloat(ang)) - # if we have a specified angular range, use that - if len(args) > 0: + # if we have a specified angular range, use it + if range is not None: angRange = np.atleast_1d(npfloat(args[0])) # divide of multiples of period @@ -549,6 +544,11 @@ def map_angle(ang, *args, **kwargs): ubi = ang > ub pass retval = ang + # shouldn't all this be equivalent to: + # retval = np.mod(ang - lb, period) + lb ???? + # note the particular case below for range (-0.5*period, +0.5*period) + # where lb would be -0.5*period. + else: retval = np.mod(ang + 0.5*period, period) - 0.5*period return retval @@ -564,8 +564,8 @@ def row_norm(a): @xf_api -def unit_vector(a): - a = np.atleast_2d(a) +def unit_vector(vec_in): + a = np.atleast_2d(vec_in) n = a.shape[1] # calculate row norms and prevent divide by zero @@ -619,9 +619,9 @@ def make_rmat_of_expmap(exp_map): @xf_api -def make_binary_rmat(n): - assert len(n) == 3, 'Axis input does not have 3 components' - return 2*np.outer(n, n) - cnst.identity_3x3 +def make_binary_rmat(axis): + assert len(axis) == 3, 'Axis input does not have 3 components' + return 2*np.outer(axis, axis) - cnst.identity_3x3 @xf_api @@ -660,14 +660,14 @@ def angles_in_range(angles, starts, stops, degrees=True): @xf_api -def validate_angle_ranges(ang_list, startAngs, stopAngs, ccw=True): +def validate_angle_ranges(ang_list, start_angs, stop_angs, ccw=True): # Prefer ravel over flatten because flatten never skips the copy ang_list = np.asarray(ang_list).ravel() - startAngs = np.asarray(startAngs).ravel() - stopAngs = np.asarray(stopAngs).ravel() + startAngs = np.asarray(start_angs).ravel() + stopAngs = np.asarray(stop_angs).ravel() - n_ranges = len(startAngs) - assert len(stopAngs) == n_ranges, \ + n_ranges = len(start_angs) + assert len(stop_angs) == n_ranges, \ "length of min and max angular limits must match!" # to avoid warnings in >=, <= later down, mark nans; From 6377afe5f0f5a05817aaedcfdcc1f3b090053bd1 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 15 Oct 2018 18:07:21 +0200 Subject: [PATCH 12/30] removed unneeded trash from definitions. --- hexrd/transforms/transforms_definitions.py | 115 +-------------------- 1 file changed, 2 insertions(+), 113 deletions(-) diff --git a/hexrd/transforms/transforms_definitions.py b/hexrd/transforms/transforms_definitions.py index 3535d9b3..30300bd5 100644 --- a/hexrd/transforms/transforms_definitions.py +++ b/hexrd/transforms/transforms_definitions.py @@ -86,14 +86,6 @@ def _signature(angs, rmat_c=None): pass - @property - def _args(self): - return ('angs',) - - @property - def _kwargs(self): - return ('beam_vec', 'eta_vec', 'chi', 'rmat_c') - class DEF_angles_to_dvec(DEF_Func): """ @@ -113,14 +105,6 @@ def _signature(angs, rmat_c=None): pass - @property - def _args(self): - return ('angs',) - - @property - def _kwargs(self): - return ('beam_vec', 'eta_vec', 'chi', 'rmat_c') - class DEF_gvec_to_xy(DEF_Func): """ @@ -192,15 +176,6 @@ def _signature(gvec_c, bmat=None): pass - @property - def _args(self): - return ('gvec_c', 'rmat_d', 'rmat_s', 'rmat_c', 'tvec_d', 'tvec_s', - 'tvec_c') - - @property - def _kwargs(self): - return ('beam_vec', 'vmat_inv', 'bmat') - class DEF_xy_to_gvec(DEF_Func): """ @@ -259,14 +234,6 @@ def _signature(xy_d, output_ref=False): pass - @property - def _args(self): - return ('xy_d', 'rmat_d', 'rmat_s', 'tvec_d', 'tvec_s', 'tvec_c') - - @property - def _kwargs(self): - return ('rmat_b', 'distortion', 'output_ref') - class DEF_solve_omega(DEF_Func): """ @@ -356,14 +323,6 @@ def _signature(gvecs, chi, rmat_c, wavelength, bmat=None, vmat_inv=None, rmat_b=None): pass - @property - def _args(self): - return ('gvecs', 'chi', 'rmat_c', 'wavelength') - - @property - def _kwargs(self): - return ('bmat', 'vmat_inv', 'rmat_b') - # ============================================================================== # UTILITY FUNCTIONS API @@ -378,14 +337,6 @@ class DEF_angular_difference(DEF_Func): def _signature(ang_list0, ang_list1, units=cnst.angular_units): pass - @property - def _args(self): - return ('ang_list0', 'ang_list1') - - @property - def _kwargs(self): - return ('units',) - class DEF_map_angle(DEF_Func): """ @@ -401,14 +352,6 @@ class DEF_map_angle(DEF_Func): def _signature(ang, range=None, units=cnst.angular_units): pass - @property - def _args(self): - return ('ang',) - - @property - def _kwargs(self): - return ('range', 'units') - class DEF_row_norm(DEF_Func): """ @@ -417,22 +360,16 @@ class DEF_row_norm(DEF_Func): def _signature(a): pass - @property - def _args(self): - return ('a') - class DEF_unit_vector(DEF_Func): """ normalize an array of row vectors (vstacked, axis=0) + + guaranteed to work for 1d and 2d arrays. """ def _signature(vec_in): pass - @property - def _args(self): - return('vec_in') - class DEF_make_sample_rmat(DEF_Func): """ @@ -446,10 +383,6 @@ class DEF_make_sample_rmat(DEF_Func): def _signature(chi, ome): pass - @property - def _args(self): - return ('chi', 'ome') - class DEF_make_rmat_of_expmap(DEF_Func): """ @@ -458,10 +391,6 @@ class DEF_make_rmat_of_expmap(DEF_Func): def _signature(exp_map): pass - @property - def _args(self): - return ('exp_map',) - class DEF_make_binary_rmat(DEF_Func): """ @@ -470,10 +399,6 @@ class DEF_make_binary_rmat(DEF_Func): def _signature(axis): pass - @property - def _args(self): - return ('axis',) - class DEF_make_beam_rmat(DEF_Func): """ @@ -484,10 +409,6 @@ class DEF_make_beam_rmat(DEF_Func): def _signature(bvec_l, evec_l): pass - @property - def _args(self): - return ('bvec_l', 'evec_l') - class DEF_angles_in_range(DEF_Func): """Determine whether angles lie in or out of specified ranges @@ -502,14 +423,6 @@ class DEF_angles_in_range(DEF_Func): def _signature(angles, starts, stops, degrees=True): pass - @property - def _args(self): - return ('angles', 'starts', 'stops') - - @property - def _kwargs(self): - return ('degrees',) - class DEF_validate_angle_ranges(DEF_Func): """ @@ -521,14 +434,6 @@ class DEF_validate_angle_ranges(DEF_Func): def _signature(ang_list, start_angs, stop_angs, ccw=True): pass - @property - def _args(self): - return ('ang_list', 'start_angs', 'stop_angs') - - @property - def _kwargs(self): - return ('ccw',) - class DEF_rotate_vecs_about_axis(DEF_Func): """ @@ -553,10 +458,6 @@ class DEF_rotate_vecs_about_axis(DEF_Func): def _signature(angle, axis, vecs): pass - @property - def _args(self): - return ('angle', 'axis', 'vecs') - class DEF_quat_product_matrix(DEF_Func): """ @@ -588,14 +489,6 @@ class DEF_quat_product_matrix(DEF_Func): def _signature(q, mult='right'): pass - @property - def _args(self): - return ('q',) - - @property - def _kwargs(self): - return ('mult',) - class DEF_quat_distance(DEF_Func): """ @@ -604,10 +497,6 @@ class DEF_quat_distance(DEF_Func): def _signature(q1, q2, qsym): pass - @property - def _args(self): - return ('q1', 'q2', 'qsym') - # ============================================================================== # Decorator to mark implementations of the API. Names must match. From 95872c742ce94c910d1da833562b66173a68d45d Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 15 Oct 2018 18:48:23 +0200 Subject: [PATCH 13/30] renamed test_angle_in_range to match function name (test_angles_in_range) --- .../tests/{test_angle_in_range.py => test_angles_in_range.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename hexrd/transforms/tests/{test_angle_in_range.py => test_angles_in_range.py} (100%) diff --git a/hexrd/transforms/tests/test_angle_in_range.py b/hexrd/transforms/tests/test_angles_in_range.py similarity index 100% rename from hexrd/transforms/tests/test_angle_in_range.py rename to hexrd/transforms/tests/test_angles_in_range.py From b5d91aa46e25cdc5f76a03ce2ae4df5e10b0e18c Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Wed, 24 Oct 2018 15:35:20 +0200 Subject: [PATCH 14/30] added numba versions of angles_to_gvec and make_beam_rmat MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Also added a workaround to the jit dispatcher not supporting ‘signature’ --- hexrd/transforms/tests/test_angles_to_gvec.py | 4 +- hexrd/transforms/tests/test_make_beam_rmat.py | 4 +- hexrd/transforms/xf_numba.py | 110 +++++++++++------- 3 files changed, 75 insertions(+), 43 deletions(-) diff --git a/hexrd/transforms/tests/test_angles_to_gvec.py b/hexrd/transforms/tests/test_angles_to_gvec.py index ba72cd41..90fda28e 100644 --- a/hexrd/transforms/tests/test_angles_to_gvec.py +++ b/hexrd/transforms/tests/test_angles_to_gvec.py @@ -5,14 +5,14 @@ from .. import angles_to_gvec as default_angles_to_gvec from ..xf_numpy import angles_to_gvec as numpy_angles_to_gvec from ..xf_capi import angles_to_gvec as capi_angles_to_gvec -#from ..xf_numba import angles_to_gvec as numba_angles_to_gvec +from ..xf_numba import angles_to_gvec as numba_angles_to_gvec import pytest all_impls = pytest.mark.parametrize('angles_to_gvec_impl, module_name', [(numpy_angles_to_gvec, 'numpy'), (capi_angles_to_gvec, 'capi'), - #(numba_angles_to_gvec, 'numba'), + (numba_angles_to_gvec, 'numba'), (default_angles_to_gvec, 'default')] ) diff --git a/hexrd/transforms/tests/test_make_beam_rmat.py b/hexrd/transforms/tests/test_make_beam_rmat.py index 95b9522a..2f83b973 100644 --- a/hexrd/transforms/tests/test_make_beam_rmat.py +++ b/hexrd/transforms/tests/test_make_beam_rmat.py @@ -5,14 +5,14 @@ from .. import make_beam_rmat as default_make_beam_rmat from ..xf_numpy import make_beam_rmat as numpy_make_beam_rmat from ..xf_capi import make_beam_rmat as capi_make_beam_rmat -#from ..xf_numba import make_beam_rmat as numba_make_beam_rmat +from ..xf_numba import make_beam_rmat as numba_make_beam_rmat import pytest all_impls = pytest.mark.parametrize('make_beam_rmat_impl, module_name', [(numpy_make_beam_rmat, 'numpy'), (capi_make_beam_rmat, 'capi'), - #(numba_angles_to_gvec, 'numba'), + (numba_make_beam_rmat, 'numba'), (default_make_beam_rmat, 'default')] ) diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index 3dbbb820..91e3c31d 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -35,9 +35,18 @@ # from hexrd import constants as cnst from .. import constants as cnst +from .transforms_definitions import xf_api, get_signature import numba +# Use the following decorator instead of numba.jit for interface functions. +# This is so we can patch certain features. +def xfapi_jit(fn): + out = numba.jit(fn) + out.__signature__ = get_signature(fn) + + return out + def _beam_to_crystal(vecs, rmat_b=None, rmat_s=None, rmat_c=None): """ @@ -81,15 +90,18 @@ def _beam_to_crystal(vecs, rmat_b=None, rmat_s=None, rmat_c=None): @numba.njit -def _angles_to_gvec_helper(angs, out): +def _angles_to_gvec_helper(angs): """ angs are vstacked [2*theta, eta, omega] - gvec_b = np.vstack([[np.cos(0.5*angs[:, 0]) * np.cos(angs[:, 1])], - [np.cos(0.5*angs[:, 0]) * np.sin(angs[:, 1])], - [np.sin(0.5*angs[:, 0])]]) + This should be equivalent to the one-liner numpy version: + out = np.vstack([[np.cos(0.5*angs[:, 0]) * np.cos(angs[:, 1])], + [np.cos(0.5*angs[:, 0]) * np.sin(angs[:, 1])], + [np.sin(0.5*angs[:, 0])]]) + + although much faster, """ - n = angs.shape[0] + out = np.empty((nvecs, 3)) for i in range(n): ca0 = np.cos(0.5*angs[i, 0]) sa0 = np.sin(0.5*angs[i, 0]) @@ -99,6 +111,8 @@ def _angles_to_gvec_helper(angs, out): out[i, 1] = ca0 * sa1 out[i, 2] = sa0 + return out + @numba.njit def _angles_to_dvec_helper(angs, out): @@ -143,35 +157,42 @@ def _rmat_s_helper(chi, omes, out): out[i, 2, 2] = cx * cw -def angles_to_gvec(angs, rmat_b=None, chi=None, rmat_c=None): - """ - Takes triplets of angles in the beam frame (2*theta, eta, omega) - to components of unit G-vectors in the LAB frame. If the omega - values are not trivial (i.e. angs[:, 2] = 0.), then the components - are in the SAMPLE frame. If the crystal rmat is specified and - is not the identity, then the components are in the CRYSTAL frame. +@xf_api +def angles_to_gvec(angs, + beam_vec=None, eta_vec=None, + chi=None, rmat_c=None): + """Note about this implementation: + This used to take rmat_b instead of the pair beam_vec, eta_vec. So it may require + some checking. """ + # apply defaults to beam_vec and eta_vec: + beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec + eta_vec = eta_vec if eta_vec is not None else cnst.beam_vec + angs = np.atleast_2d(angs) nvecs, dim = angs.shape # make vectors in beam frame - gvec_b = np.empty((nvecs, 3)) - _angles_to_gvec_helper(angs, gvec_b) + gvec_b = _angles_to_gvec_helper(angs) + + # for NUMBA case, need chi as a flaot if chi is None: chi = 0. - # calcuate rmat_s + # calculate rmat_s if dim > 2: rmat_s = np.empty((nvecs, 3, 3)) _rmat_s_helper(chi, angs[:, 2], rmat_s) else: # no omegas given rmat_s = np.empty((1, 3, 3)) _rmat_s_helper(chi, [0., ], rmat_s) - return _beam_to_crystal( - gvec_b, rmat_b=rmat_b, rmat_s=rmat_s, rmat_c=rmat_c - ) + + rmat_b = make_beam_rmat(beam_vec, eta_vec) + out = _beam_to_crystal(gvec_b, + rmat_b=rmat_b, rmat_s=rmat_s, rmat_c=rmat_c) + return out def angles_to_dvec(angs, rmat_b=None, chi=None, rmat_c=None): @@ -336,44 +357,55 @@ def make_rmat_of_expmap(exp_map): return np.squeeze(rmats) -@numba.njit -def _make_beam_rmat(bhat_l, ehat_l, out): +@xf_api +@xfapi_jit +def make_beam_rmat(bvec_l, evec_l): # bhat_l and ehat_l CANNOT have 0 magnitude! # must catch this case as well as colinear bhat_l/ehat_l elsewhere... - bHat_mag = np.sqrt(bhat_l[0]**2 + bhat_l[1]**2 + bhat_l[2]**2) + + bvec_mag = np.sqrt(bvec_l[0]**2 + bvec_l[1]**2 + bvec_l[2]**2) + + if bvec_mag < cnst.sqrt_epsf: + #can numba raise? + # raise RuntimeError("beam_vec MUST NOT be ZERO!") + pass # assign Ze as -bhat_l for i in range(3): - out[i, 2] = -bhat_l[i] / bHat_mag + out[i, 2] = -bvec_l[i] / bvec_mag + Ze0 = -bvec_l[0] / bvec_mag + Ze1 = -bvec_l[1] / bvec_mag + Ze2 = -bvec_l[2] / bvec_mag # find Ye as Ze ^ ehat_l - Ye0 = out[1, 2]*ehat_l[2] - ehat_l[1]*out[2, 2] - Ye1 = out[2, 2]*ehat_l[0] - ehat_l[2]*out[0, 2] - Ye2 = out[0, 2]*ehat_l[1] - ehat_l[0]*out[1, 2] + Ye0 = Ze1*evec_l[2] - evec_l[1]*Ze2 + Ye1 = Ze2*evec_l[0] - evec_l[2]*Ze0 + Ye2 = Ze0*evec_l[1] - evec_l[0]*Ze1 Ye_mag = np.sqrt(Ye0**2 + Ye1**2 + Ye2**2) + if Ye_mag < cnst.sqrt_epsf: + # raise RuntimeError("beam_vec and eta_vec MUST NOT be colinear!") + pass + + out = np.empty((3,3)) # numba can now allocate + Ye0 /= Ye_mag + Ye1 /= Ye_mag + Ye2 /= Ye_mag + + # find Xe as Ye ^ Ze + Xe0 = Ye1*Ze2 - Ze1*Ye2 + Xe1 = Ye2*Ze0 - Ze2*Ye0 + Xe2 = Ye0*Ze1 - Ze0*Ye1 + out[0, 1] = Ye0 / Ye_mag out[1, 1] = Ye1 / Ye_mag out[2, 1] = Ye2 / Ye_mag - # find Xe as Ye ^ Ze out[0, 0] = out[1, 1]*out[2, 2] - out[1, 2]*out[2, 1] out[1, 0] = out[2, 1]*out[0, 2] - out[2, 2]*out[0, 1] out[2, 0] = out[0, 1]*out[1, 2] - out[0, 2]*out[1, 1] + return out -def make_beam_rmat(bhat_l, ehat_l): - """ - make eta basis COB matrix with beam antiparallel with Z - - takes components from BEAM frame to LAB - - FIXME: NO EXCEPTION HANDLING FOR COLINEAR ARGS IN NUMBA VERSION! - - ???: put checks for non-zero magnitudes and non-colinearity in wrapper? - """ - result = np.empty((3, 3)) - _make_beam_rmat(bhat_l.reshape(3), ehat_l.reshape(3), result) - return result From 43103f13e8e8555596c615551ce77eff08e43981 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Wed, 24 Oct 2018 19:12:12 +0200 Subject: [PATCH 15/30] added angles_to_dvec to numba implemented includes some clean up of numba functions. --- hexrd/transforms/tests/test_angles_to_dvec.py | 4 +- hexrd/transforms/xf_numba.py | 185 +++++++++++------- 2 files changed, 117 insertions(+), 72 deletions(-) diff --git a/hexrd/transforms/tests/test_angles_to_dvec.py b/hexrd/transforms/tests/test_angles_to_dvec.py index f4d32432..293f64e0 100644 --- a/hexrd/transforms/tests/test_angles_to_dvec.py +++ b/hexrd/transforms/tests/test_angles_to_dvec.py @@ -5,14 +5,14 @@ from .. import angles_to_dvec as default_angles_to_dvec from ..xf_numpy import angles_to_dvec as numpy_angles_to_dvec from ..xf_capi import angles_to_dvec as capi_angles_to_dvec -#from ..xf_numba import angles_to_dvec as numba_angles_to_dvec +from ..xf_numba import angles_to_dvec as numba_angles_to_dvec import pytest all_impls = pytest.mark.parametrize('angles_to_dvec_impl, module_name', [(numpy_angles_to_dvec, 'numpy'), (capi_angles_to_dvec, 'capi'), - #(numba_angles_to_gvec, 'numba'), + (numba_angles_to_dvec, 'numba'), (default_angles_to_dvec, 'default')] ) diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index 91e3c31d..f815c94a 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -90,19 +90,21 @@ def _beam_to_crystal(vecs, rmat_b=None, rmat_s=None, rmat_c=None): @numba.njit -def _angles_to_gvec_helper(angs): +def _angles_to_gvec_helper(angs, out=None): """ - angs are vstacked [2*theta, eta, omega] + angs are vstacked [2*theta, eta, omega], although omega is optional This should be equivalent to the one-liner numpy version: out = np.vstack([[np.cos(0.5*angs[:, 0]) * np.cos(angs[:, 1])], [np.cos(0.5*angs[:, 0]) * np.sin(angs[:, 1])], [np.sin(0.5*angs[:, 0])]]) - although much faster, + although much faster """ - out = np.empty((nvecs, 3)) - for i in range(n): + _, dim = angs.shape + out = out if out is not None else np.empty((dim, 3), dtype=angs.dtype) + + for i in range(len(angs)): ca0 = np.cos(0.5*angs[i, 0]) sa0 = np.sin(0.5*angs[i, 0]) ca1 = np.cos(angs[i, 1]) @@ -115,16 +117,20 @@ def _angles_to_gvec_helper(angs): @numba.njit -def _angles_to_dvec_helper(angs, out): +def _angles_to_dvec_helper(angs, out=None): """ - angs are vstacked [2*theta, eta, omega] + angs are vstacked [2*theta, eta, omega], although omega is optional + + This shoud be equivalent to the one-liner numpy version: + out = np.vstack([[np.sin(angs[:, 0]) * np.cos(angs[:, 1])], + [np.sin(angs[:, 0]) * np.sin(angs[:, 1])], + [-np.cos(angs[:, 0])]]) - dvec_b = np.vstack([[np.sin(angs[:, 0]) * np.cos(angs[:, 1])], - [np.sin(angs[:, 0]) * np.sin(angs[:, 1])], - [-np.cos(angs[:, 0])]]) + although much faster """ - n = angs.shape[0] - for i in range(n): + _, dim = angs.shape + out = out if out is not None else np.empty((dim, 3), dtype=angs.dtype) + for i in range(len(angs)): ca0 = np.cos(angs[i, 0]) sa0 = np.sin(angs[i, 0]) ca1 = np.cos(angs[i, 1]) @@ -133,28 +139,71 @@ def _angles_to_dvec_helper(angs, out): out[i, 1] = sa0 * sa1 out[i, 2] = -ca0 + return out @numba.njit -def _rmat_s_helper(chi, omes, out): +def _rmat_s_helper(chi=None, omes=None, out=None): """ simple utility for calculating sample rotation matrices based on standard definition for HEDM + + chi is a single value, 0.0 by default + omes is either a 1d array or None. + If None the code should be equivalent to a single ome of value 0.0 + + out is a preallocated output array. No check is done about it having the + proper size. If None a new array will be allocated. The expected size + of the array is as many 3x3 matrices as omes (n, 3, 3). """ - n = len(omes) - cx = np.cos(chi) - sx = np.sin(chi) - for i in range(n): - cw = np.cos(omes[i]) - sw = np.sin(omes[i]) - out[i, 0, 0] = cw - out[i, 0, 1] = 0. - out[i, 0, 2] = sw - out[i, 1, 0] = sx * sw - out[i, 1, 1] = cx - out[i, 1, 2] = -sx * cw - out[i, 2, 0] = -cx * sw - out[i, 2, 1] = sx - out[i, 2, 2] = cx * cw + if chi is not None: + cx = np.cos(chi) + sx = np.sin(chi) + else: + cx = 1.0 + sx = 0.0 + + if omes is not None: + # omes is an array (vector): output is as many rotation matrices as omes entries. + n = len(omes) + out = out if out is not None else np.empty((n,3,3), dtype=omes.dtype) + + if chi is not None: + # ome is array and chi is a value... compute output + cx = np.cos(chi) + sx = np.sin(chi) + for i in range(n): + cw = np.cos(omes[i]) + sw = np.sin(omes[i]) + out[i, 0, 0] = cw; out[i, 0, 1] = 0.; out[i, 0, 2] = sw + out[i, 1, 0] = sx*sw; out[i, 1, 1] = cx; out[i, 1, 2] = -sx*cw + out[i, 2, 0] = -cx*sw; out[i, 2, 1] = sx; out[i, 2, 2] = cx*cw + else: + # omes is array and chi is None -> equivalent to chi=0.0, but shortcut computations. + # cx IS 1.0, sx IS 0.0 + for i in range(n): + cw = np.cos(omes[i]) + sw = np.sin(omes[i]) + out[i, 0, 0] = cw; out[i, 0, 1] = 0.; out[i, 0, 2] = sw + out[i, 1, 0] = 0.; out[i, 1, 1] = 1.; out[i, 1, 2] = 0. + out[i, 2, 0] = -sw; out[i, 2, 1] = 0.; out[i, 2, 2] = cw + else: + # omes is None, results should be equivalent to an array with a single element 0.0 + out = out if out is not None else np.empty((1, 3, 3)) + if chi is not None: + # ome is 0.0. cw is 1.0 and sw is 0.0 + cx = np.cos(chi) + sx = np.sin(chi) + out[0, 0, 0] = 1.; out[0, 0, 1] = 0.; out[0, 0, 2] = 0. + out[0, 1, 0] = 0.; out[0, 1, 1] = cx; out[0, 1, 2] = -sx + out[0, 2, 0] = 0.; out[0, 2, 1] = sx; out[0, 2, 2] = cx + else: + # both omes and chi are None... return a single identity matrix. + out[0, 0, 0] = 1.; out[0, 0, 1] = 0.; out[0, 0, 2] = 0. + out[0, 1, 0] = 0.; out[0, 1, 1] = 1.; out[0, 1, 2] = 0. + out[0, 2, 0] = 0.; out[0, 2, 1] = 0.; out[0, 2, 2] = 1. + + + return out @xf_api @@ -165,65 +214,61 @@ def angles_to_gvec(angs, This used to take rmat_b instead of the pair beam_vec, eta_vec. So it may require some checking. """ - # apply defaults to beam_vec and eta_vec: - beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec - eta_vec = eta_vec if eta_vec is not None else cnst.beam_vec - angs = np.atleast_2d(angs) nvecs, dim = angs.shape # make vectors in beam frame - gvec_b = _angles_to_gvec_helper(angs) - - + gvec_b = _angles_to_gvec_helper(angs[:,0:2]) - # for NUMBA case, need chi as a flaot - if chi is None: - chi = 0. - - # calculate rmat_s - if dim > 2: - rmat_s = np.empty((nvecs, 3, 3)) - _rmat_s_helper(chi, angs[:, 2], rmat_s) - else: # no omegas given - rmat_s = np.empty((1, 3, 3)) - _rmat_s_helper(chi, [0., ], rmat_s) + # _rmat_s_helper could return None to mean "Identity" when chi and ome are None. + omes = angs[:, 2] if dim>2 else None + if chi is not None or omes is not None: + rmat_s = _rmat_s_helper(chi=chi, omes=omes) + else: + rmat_s = None + # apply defaults to beam_vec and eta_vec. + # TODO: use a default rmat when beam_vec and eta_vec are None so computations + # can be avoided? + beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec + eta_vec = eta_vec if eta_vec is not None else cnst.beam_vec rmat_b = make_beam_rmat(beam_vec, eta_vec) + out = _beam_to_crystal(gvec_b, rmat_b=rmat_b, rmat_s=rmat_s, rmat_c=rmat_c) return out -def angles_to_dvec(angs, rmat_b=None, chi=None, rmat_c=None): - """ - Takes triplets of angles in the beam frame (2*theta, eta, omega) - to components of unit diffraction vectors in the LAB frame. If the - omega values are not trivial (i.e. angs[:, 2] = 0.), then the - components are in the SAMPLE frame. If the crystal rmat is specified - and is not the identity, then the components are in the CRYSTAL frame. +@xf_api +def angles_to_dvec(angs, + beam_vec=None, eta_vec=None, + chi=None, rmat_c=None): + """Note about this implementation: + This used to take rmat_b instead of the pair beam_vec, eta_vec. So it may require + some checking. """ angs = np.atleast_2d(angs) nvecs, dim = angs.shape # make vectors in beam frame - dvec_b = np.empty((nvecs, 3)) - _angles_to_dvec_helper(angs, dvec_b) - - # for NUMBA case, need chi as a flaot - if chi is None: - chi = 0. - - # calcuate rmat_s - if dim > 2: - rmat_s = np.empty((nvecs, 3, 3)) - _rmat_s_helper(chi, angs[:, 2], rmat_s) - else: # no omegas given - rmat_s = np.empty((1, 3, 3)) - _rmat_s_helper(chi, [0., ], rmat_s) - return _beam_to_crystal( - dvec_b, rmat_b=rmat_b, rmat_s=rmat_s, rmat_c=rmat_c - ) + dvec_b = _angles_to_dvec_helper(angs[:,0:2]) + + # calculate rmat_s + omes = angs[:, 2] if dim>2 else None + if chi is not None or omes is not None: + rmat_s = _rmat_s_helper(chi=chi, omes=omes) + else: + rmat_s = None + + # apply defaults to beam_vec and eta_vec. + # TODO: use a default rmat when beam_vec and eta_vec are None so computations + # can be avoided? + beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec + eta_vec = eta_vec if eta_vec is not None else cnst.beam_vec + rmat_b = make_beam_rmat(beam_vec, eta_vec) + + return _beam_to_crystal(dvec_b, + rmat_b=rmat_b, rmat_s=rmat_s, rmat_c=rmat_c) @numba.njit From ac18aae7d9f398c3926b93cb269aad6f01fb392c Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Wed, 24 Oct 2018 19:28:25 +0200 Subject: [PATCH 16/30] added numba row_norm --- hexrd/transforms/tests/test_row_norm.py | 4 ++-- hexrd/transforms/xf_numba.py | 24 +++++++++++++++++------- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/hexrd/transforms/tests/test_row_norm.py b/hexrd/transforms/tests/test_row_norm.py index 55e440b5..5675a05f 100644 --- a/hexrd/transforms/tests/test_row_norm.py +++ b/hexrd/transforms/tests/test_row_norm.py @@ -5,14 +5,14 @@ from .. import row_norm as default_row_norm from ..xf_numpy import row_norm as numpy_row_norm #from ..xf_capi import row_norm as capi_row_norm -#from ..xf_numba import row_norm as numba_row_norm +from ..xf_numba import row_norm as numba_row_norm import pytest all_impls = pytest.mark.parametrize('row_norm_impl, module_name', [(numpy_row_norm, 'numpy'), #(capi_row_norm, 'capi'), - #(numba_angles_to_gvec, 'numba'), + (numba_row_norm, 'numba'), (default_row_norm, 'default')] ) diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index f815c94a..c8d6a747 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -271,14 +271,19 @@ def angles_to_dvec(angs, rmat_b=rmat_b, rmat_s=rmat_s, rmat_c=rmat_c) +# this could be a gufunc... (n)->() @numba.njit -def _row_norm(a, b): +def _row_norm(a, out=None): n, dim = a.shape + out = out if out is not None else np.empty(n, dtype=a.dtype) for i in range(n): nrm = 0.0 for j in range(dim): - nrm += a[i, j]*a[i, j] - b[i] = np.sqrt(nrm) + x = a[i, j] + nrm += x*x + out[i] = np.sqrt(nrm) + + return out @numba.njit @@ -314,14 +319,19 @@ def _unit_vector_multi(a, b): b[i, j] = a[i, j] +@xf_api def row_norm(a): """ - retrun row-wise norms for a list of vectors + return row-wise norms for a list of vectors """ + # TODO: leave this to a PRECONDITION in the xf_api? + if len(a.shape)>2: + raise RuntimeError( + "incorrect shape: arg must be 1-d or 2-d, yours is %d" + % (len(a.shape))) + a = np.atleast_2d(a) - result = np.empty(len(a)) - _row_norm(a, result) - return result + return _row_norm(a, result) def unit_vector(a): From 4898e99b1b6963dff886a16fb2d53e7417f7fda2 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Wed, 24 Oct 2018 20:11:09 +0200 Subject: [PATCH 17/30] added numba versions of unit_vector and make_rmat_of_expmap --- .../tests/test_make_rmat_of_expmap.py | 4 +- hexrd/transforms/tests/test_unit_vector.py | 4 +- hexrd/transforms/xf_numba.py | 83 ++++++++++--------- 3 files changed, 48 insertions(+), 43 deletions(-) diff --git a/hexrd/transforms/tests/test_make_rmat_of_expmap.py b/hexrd/transforms/tests/test_make_rmat_of_expmap.py index ac6b4074..b33bacaa 100644 --- a/hexrd/transforms/tests/test_make_rmat_of_expmap.py +++ b/hexrd/transforms/tests/test_make_rmat_of_expmap.py @@ -5,14 +5,14 @@ from .. import make_rmat_of_expmap as default_make_rmat_of_expmap from ..xf_numpy import make_rmat_of_expmap as numpy_make_rmat_of_expmap from ..xf_capi import make_rmat_of_expmap as capi_make_rmat_of_expmap -#from ..xf_numba import make_rmat_of_expmap as numba_make_rmat_of_expmap +from ..xf_numba import make_rmat_of_expmap as numba_make_rmat_of_expmap import pytest all_impls = pytest.mark.parametrize('make_rmat_of_expmap_impl, module_name', [(numpy_make_rmat_of_expmap, 'numpy'), (capi_make_rmat_of_expmap, 'capi'), - #(numba_angles_to_gvec, 'numba'), + (numba_make_rmat_of_expmap, 'numba'), (default_make_rmat_of_expmap, 'default')] ) diff --git a/hexrd/transforms/tests/test_unit_vector.py b/hexrd/transforms/tests/test_unit_vector.py index bf140840..47d85617 100644 --- a/hexrd/transforms/tests/test_unit_vector.py +++ b/hexrd/transforms/tests/test_unit_vector.py @@ -5,14 +5,14 @@ from .. import unit_vector as default_unit_vector from ..xf_numpy import unit_vector as numpy_unit_vector from ..xf_capi import unit_vector as capi_unit_vector -#from ..xf_numba import unit_vector as numba_unit_vector +from ..xf_numba import unit_vector as numba_unit_vector import pytest all_impls = pytest.mark.parametrize('unit_vector_impl, module_name', [(numpy_unit_vector, 'numpy'), (capi_unit_vector, 'capi'), - #(numba_angles_to_gvec, 'numba'), + (numba_unit_vector, 'numba'), (default_unit_vector, 'default')] ) diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index c8d6a747..70598876 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -286,8 +286,11 @@ def _row_norm(a, out=None): return out +# this and _unit_vector_single would be better as a gufunc. @numba.njit -def _unit_vector_single(a, b): +def _unit_vector_single(a, out=None): + out = out if out is not None else np.empty_like(a) + n = len(a) nrm = 0.0 for i in range(n): @@ -296,14 +299,17 @@ def _unit_vector_single(a, b): # prevent divide by zero if nrm > cnst.epsf: for i in range(n): - b[i] = a[i] / nrm + out[i] = a[i] / nrm else: for i in range(n): - b[i] = a[i] + out[i] = a[i] + return out @numba.njit -def _unit_vector_multi(a, b): +def _unit_vector_multi(a, out=None): + out = out if out is not None else np.empty_like(a) + n, dim = a.shape for i in range(n): nrm = 0.0 @@ -313,11 +319,12 @@ def _unit_vector_multi(a, b): # prevent divide by zero if nrm > cnst.epsf: for i in range(n): - b[i, j] = a[i, j] / nrm + out[i, j] = a[i, j] / nrm else: for i in range(n): - b[i, j] = a[i, j] + out[i, j] = a[i, j] + return out @xf_api def row_norm(a): @@ -331,28 +338,27 @@ def row_norm(a): % (len(a.shape))) a = np.atleast_2d(a) - return _row_norm(a, result) + return _row_norm(a) - -def unit_vector(a): +@xf_api +def unit_vector(vec_in): """ normalize array of column vectors (hstacked, axis = 0) """ - result = np.empty_like(a) - if a.ndim == 1: - _unit_vector_single(a, result) - elif a.ndim == 2: - _unit_vector_multi(a, result) + if vec_in.ndim == 1: + out = _unit_vector_single(vec_in) + elif vec_in.ndim == 2: + out = _unit_vector_multi(vec_in) else: raise ValueError( "incorrect arg shape; must be 1-d or 2-d, yours is %d-d" - % (a.ndim) + % (vec_in.ndim) ) - return result + return out @numba.njit -def _make_rmat_of_expmap(x, z): +def _make_rmat_of_expmap(x, out=None): """ TODO: @@ -362,31 +368,31 @@ def _make_rmat_of_expmap(x, z): for the phi = 0 cases, and deal with it later; or 2) catch phi = 0 cases inside the loop and just return squeezed answer """ - n = x.shape[0] + n = len(x) + out = out if out is not None else np.empty((n,3,3), dtype=x.dtype) for i in range(n): phi = np.sqrt(x[i, 0]*x[i, 0] + x[i, 1]*x[i, 1] + x[i, 2]*x[i, 2]) if phi <= cnst.sqrt_epsf: - z[i, 0, 0] = 1. - z[i, 0, 1] = 0. - z[i, 0, 2] = 0. - z[i, 1, 0] = 0. - z[i, 1, 1] = 1. - z[i, 1, 2] = 0. - z[i, 2, 0] = 0. - z[i, 2, 1] = 0. - z[i, 2, 2] = 1. + out[i, 0, 0] = 1.; out[i, 0, 1] = 0.; out[i, 0, 2] = 0. + out[i, 1, 0] = 0.; out[i, 1, 1] = 1.; out[i, 1, 2] = 0. + out[i, 2, 0] = 0.; out[i, 2, 1] = 0.; out[i, 2, 2] = 1. else: f1 = np.sin(phi)/phi f2 = (1. - np.cos(phi)) / (phi*phi) - z[i, 0, 0] = 1. - f2*(x[i, 2]*x[i, 2] + x[i, 1]*x[i, 1]) - z[i, 0, 1] = f2*x[i, 1]*x[i, 0] - f1*x[i, 2] - z[i, 0, 2] = f1*x[i, 1] + f2*x[i, 2]*x[i, 0] - z[i, 1, 0] = f1*x[i, 2] + f2*x[i, 1]*x[i, 0] - z[i, 1, 1] = 1. - f2*(x[i, 2]*x[i, 2] + x[i, 0]*x[i, 0]) - z[i, 1, 2] = f2*x[i, 2]*x[i, 1] - f1*x[i, 0] - z[i, 2, 0] = f2*x[i, 2]*x[i, 0] - f1*x[i, 1] - z[i, 2, 1] = f1*x[i, 0] + f2*x[i, 2]*x[i, 1] - z[i, 2, 2] = 1. - f2*(x[i, 1]*x[i, 1] + x[i, 0]*x[i, 0]) + + out[i, 0, 0] = 1. - f2*(x[i, 2]*x[i, 2] + x[i, 1]*x[i, 1]) + out[i, 0, 1] = f2*x[i, 1]*x[i, 0] - f1*x[i, 2] + out[i, 0, 2] = f1*x[i, 1] + f2*x[i, 2]*x[i, 0] + + out[i, 1, 0] = f1*x[i, 2] + f2*x[i, 1]*x[i, 0] + out[i, 1, 1] = 1. - f2*(x[i, 2]*x[i, 2] + x[i, 0]*x[i, 0]) + out[i, 1, 2] = f2*x[i, 2]*x[i, 1] - f1*x[i, 0] + + out[i, 2, 0] = f2*x[i, 2]*x[i, 0] - f1*x[i, 1] + out[i, 2, 1] = f1*x[i, 0] + f2*x[i, 2]*x[i, 1] + out[i, 2, 2] = 1. - f2*(x[i, 1]*x[i, 1] + x[i, 0]*x[i, 0]) + + return out """ @@ -404,11 +410,10 @@ def make_rmat_of_expmap(exp_map): return rmats """ - +@xf_api def make_rmat_of_expmap(exp_map): exp_map = np.atleast_2d(exp_map) - rmats = np.empty((len(exp_map), 3, 3)) - _make_rmat_of_expmap(exp_map, rmats) + rmats = _make_rmat_of_expmap(exp_map) return np.squeeze(rmats) From 16fe1a2395c06714dc36ed336f5b7a0d234e0746 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Wed, 24 Oct 2018 20:21:04 +0200 Subject: [PATCH 18/30] added a comment about current status of numba implementation --- hexrd/transforms/xf_numba.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index 70598876..4d6975c5 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -25,8 +25,22 @@ # the Free Software Foundation, Inc., 59 Temple Place, Suite 330, # Boston, MA 02111-1307 USA or visit . # ============================================================================= - # ??? do we want to set np.seterr(invalid='ignore') to avoid nan warnings? + +# -*- coding: utf-8 -*- +"""Tranforms module implementation using numba. + +Currently, this implementation contains code for the following functions: + +- angles_to_gvec +- angles_to_dvec + +- row_norm +- unit_vector +- make_rmat_of_expmap +- make_beam_rmap +""" + from __future__ import absolute_import import numpy as np From cf92c261d5942ffdd807b7e2d20c2ac59f65421d Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Thu, 25 Oct 2018 14:59:41 +0200 Subject: [PATCH 19/30] test in make_beam_rmat that checks an identity results from defaults. For the default beam vector and the default eta vector an identity matrix should return. This is kind of assumed in other parts of the code. --- hexrd/constants.py | 3 ++ hexrd/transforms/tests/test_make_beam_rmat.py | 22 ++++++++++---- hexrd/transforms/xf_capi.py | 5 ++-- hexrd/transforms/xf_numba.py | 30 +++++++++---------- 4 files changed, 38 insertions(+), 22 deletions(-) diff --git a/hexrd/constants.py b/hexrd/constants.py index c9f36aba..41e2dbd2 100644 --- a/hexrd/constants.py +++ b/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 diff --git a/hexrd/transforms/tests/test_make_beam_rmat.py b/hexrd/transforms/tests/test_make_beam_rmat.py index 2f83b973..aa4b348f 100644 --- a/hexrd/transforms/tests/test_make_beam_rmat.py +++ b/hexrd/transforms/tests/test_make_beam_rmat.py @@ -7,6 +7,10 @@ from ..xf_capi import make_beam_rmat as capi_make_beam_rmat from ..xf_numba import make_beam_rmat as numba_make_beam_rmat +from ... import constants as cnst + +from numpy.testing import assert_allclose + import pytest all_impls = pytest.mark.parametrize('make_beam_rmat_impl, module_name', @@ -18,9 +22,17 @@ @all_impls -def test_sample1(make_beam_rmat_impl, module_name): - pass +def test_reference_beam_mat(make_beam_rmat_impl, module_name): + """Building from the standard beam_vec and eta_vec should + yield an identity matrix. -@all_impls -def test_sample2(make_beam_rmat_impl, module_name): - pass + This is somehow assumed in other parts of the code where using the default + cnst.beam_vec and cnst.eta_vec implies an identity beam rotation matrix that + is ellided in operations""" + + beam_vec = cnst.beam_vec + eta_vec = cnst.eta_vec + + rmat = make_beam_rmat_impl(cnst.beam_vec, cnst.eta_vec) + + assert_allclose(rmat, cnst.identity_3x3) diff --git a/hexrd/transforms/xf_capi.py b/hexrd/transforms/xf_capi.py index b08e7100..e73fb3a0 100644 --- a/hexrd/transforms/xf_capi.py +++ b/hexrd/transforms/xf_capi.py @@ -27,6 +27,7 @@ from hexrd import constants as cnst from .. import _transforms_CAPI from .transforms_definitions import xf_api +import numpy as np @xf_api def angles_to_gvec( @@ -209,8 +210,8 @@ def make_binary_rmat(axis): @xf_api def make_beam_rmat(bvec_l, evec_l): - arg1 = np.ascontiguousarray(bHat_l.flatten()) - arg2 = np.ascontiguousarray(eHat_l.flatten()) + arg1 = np.ascontiguousarray(bvec_l.flatten()) + arg2 = np.ascontiguousarray(evec_l.flatten()) return _transforms_CAPI.makeEtaFrameRotMat(arg1, arg2) diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index 4d6975c5..063dc9fb 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -38,7 +38,7 @@ - row_norm - unit_vector - make_rmat_of_expmap -- make_beam_rmap +- make_beam_rmat """ from __future__ import absolute_import @@ -434,19 +434,15 @@ def make_rmat_of_expmap(exp_map): @xf_api @xfapi_jit def make_beam_rmat(bvec_l, evec_l): - # bhat_l and ehat_l CANNOT have 0 magnitude! + # bvec_l and evec_l CANNOT have 0 magnitude! # must catch this case as well as colinear bhat_l/ehat_l elsewhere... - bvec_mag = np.sqrt(bvec_l[0]**2 + bvec_l[1]**2 + bvec_l[2]**2) if bvec_mag < cnst.sqrt_epsf: - #can numba raise? - # raise RuntimeError("beam_vec MUST NOT be ZERO!") + raise RuntimeError("beam_vec MUST NOT be ZERO!") pass # assign Ze as -bhat_l - for i in range(3): - out[i, 2] = -bvec_l[i] / bvec_mag Ze0 = -bvec_l[0] / bvec_mag Ze1 = -bvec_l[1] / bvec_mag Ze2 = -bvec_l[2] / bvec_mag @@ -458,10 +454,10 @@ def make_beam_rmat(bvec_l, evec_l): Ye_mag = np.sqrt(Ye0**2 + Ye1**2 + Ye2**2) if Ye_mag < cnst.sqrt_epsf: - # raise RuntimeError("beam_vec and eta_vec MUST NOT be colinear!") + raise RuntimeError("beam_vec and eta_vec MUST NOT be colinear!") pass - out = np.empty((3,3)) # numba can now allocate + out = np.empty((3,3), dtype=bvec_l.dtype) Ye0 /= Ye_mag Ye1 /= Ye_mag Ye2 /= Ye_mag @@ -472,13 +468,17 @@ def make_beam_rmat(bvec_l, evec_l): Xe2 = Ye0*Ze1 - Ze0*Ye1 - out[0, 1] = Ye0 / Ye_mag - out[1, 1] = Ye1 / Ye_mag - out[2, 1] = Ye2 / Ye_mag + out[0, 0] = Xe0 + out[0, 1] = Ye0 + out[0, 2] = Ze0 + + out[1, 0] = Xe1 + out[1, 1] = Ye1 + out[1, 2] = Ze1 - out[0, 0] = out[1, 1]*out[2, 2] - out[1, 2]*out[2, 1] - out[1, 0] = out[2, 1]*out[0, 2] - out[2, 2]*out[0, 1] - out[2, 0] = out[0, 1]*out[1, 2] - out[0, 2]*out[1, 1] + out[2, 0] = Xe2 + out[2, 1] = Ye2 + out[2, 2] = Ze2 return out From bf326451a0e1547e199e5fd4aa050118fbb040c0 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Thu, 25 Oct 2018 16:49:59 +0200 Subject: [PATCH 20/30] tests for make_beam_rmat error conditions. Fixed api version so it generates the exceptions as the other versions do. Made the error messages a bit more coherent. --- hexrd/transforms/src/transforms_CAPI.c | 27 +++++++++++++-- hexrd/transforms/src/transforms_CFUNC.c | 34 +++++++++++++++---- hexrd/transforms/src/transforms_CFUNC.h | 9 +++-- hexrd/transforms/tests/test_make_beam_rmat.py | 21 +++++++++--- hexrd/transforms/xf_numba.py | 4 +-- hexrd/transforms/xf_numpy.py | 2 +- 6 files changed, 79 insertions(+), 18 deletions(-) diff --git a/hexrd/transforms/src/transforms_CAPI.c b/hexrd/transforms/src/transforms_CAPI.c index d3ce2934..4ac6d219 100644 --- a/hexrd/transforms/src/transforms_CAPI.c +++ b/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) diff --git a/hexrd/transforms/src/transforms_CFUNC.c b/hexrd/transforms/src/transforms_CFUNC.c index 7f15f355..b277d18c 100644 --- a/hexrd/transforms/src/transforms_CFUNC.c +++ b/hexrd/transforms/src/transforms_CFUNC.c @@ -562,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; @@ -576,10 +579,12 @@ void unitRowVector_cfunc(int n, double * cIn, double * cOut) for (j=0; j Date: Thu, 25 Oct 2018 17:51:52 +0200 Subject: [PATCH 21/30] orthonormality tests for make_beam_rmat --- hexrd/transforms/tests/test_make_beam_rmat.py | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/hexrd/transforms/tests/test_make_beam_rmat.py b/hexrd/transforms/tests/test_make_beam_rmat.py index d2f76854..3eb6c92b 100644 --- a/hexrd/transforms/tests/test_make_beam_rmat.py +++ b/hexrd/transforms/tests/test_make_beam_rmat.py @@ -49,3 +49,50 @@ def test_zero_beam_vec(make_beam_rmat_impl, module_name): def test_colinear_beam_eta_vec(make_beam_rmat_impl, module_name): with pytest.raises(RuntimeError): make_beam_rmat_impl(cnst.beam_vec, cnst.beam_vec) + + +@all_impls +def test_orthonormal_1(make_beam_rmat_impl, module_name): + beam_vec = np.array([1.0, 2.0, 3.0]) + other_vec = np.array([5.0, 2.0, 1.0]) + + # force both inputs to be orthogonal and normalized + eta_vec = np.cross(beam_vec, other_vec) + + beam_vec /= np.linalg.norm(beam_vec) + eta_vec /= np.linalg.norm(eta_vec) + + rmat = make_beam_rmat_impl(beam_vec, eta_vec) + + # dot(A, A.T) == Identity seems a good orthonormality check + # Note: atol needed as rtol is not useful for '0.' entries. + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=1e+10) + + +@all_impls +def test_orthonormal_2(make_beam_rmat_impl, module_name): + # same as above although the inputs are not normalized + beam_vec = np.array([1.0, 2.0, 3.0]) + other_vec = np.array([5.0, 2.0, 1.0]) + + # force both inputs to be orthogonal + eta_vec = np.cross(beam_vec, other_vec) + + rmat = make_beam_rmat_impl(beam_vec, eta_vec) + + # dot(A, A.T) == Identity seems a good orthonormality check + # Note: atol needed as rtol is not useful for '0.' entries. + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=1e+10) + + +@all_impls +def test_orthonormal_3(make_beam_rmat_impl, module_name): + # same as above although the inputs are neither normalized nor orthogonal + beam_vec = np.array([1.0, 2.0, 3.0]) + other_vec = np.array([5.0, 2.0, 1.0]) + + rmat = make_beam_rmat_impl(beam_vec, other_vec) + + # dot(A, A.T) == Identity seems a good orthonormality check + # Note: atol needed as rtol is not useful for '0.' entries. + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=1e+10) From dd03fe6ea90548c9abc6be12c56070ab23157925 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Fri, 26 Oct 2018 11:57:47 +0200 Subject: [PATCH 22/30] added strided argument tests for make_beam_rmat also fixed type with atol (1e+10 was meant to be 1e-10) --- hexrd/transforms/tests/test_make_beam_rmat.py | 64 +++++++++++++++++-- 1 file changed, 59 insertions(+), 5 deletions(-) diff --git a/hexrd/transforms/tests/test_make_beam_rmat.py b/hexrd/transforms/tests/test_make_beam_rmat.py index 3eb6c92b..430fc590 100644 --- a/hexrd/transforms/tests/test_make_beam_rmat.py +++ b/hexrd/transforms/tests/test_make_beam_rmat.py @@ -14,6 +14,8 @@ import pytest +ATOL_IDENTITY = 1e-10 + all_impls = pytest.mark.parametrize('make_beam_rmat_impl, module_name', [(numpy_make_beam_rmat, 'numpy'), (capi_make_beam_rmat, 'capi'), @@ -21,9 +23,12 @@ (default_make_beam_rmat, 'default')] ) +# ------------------------------------------------------------------------------ + +# Test reference frame result @all_impls -def test_reference_beam_mat(make_beam_rmat_impl, module_name): +def test_reference_beam_rmat(make_beam_rmat_impl, module_name): """Building from the standard beam_vec and eta_vec should yield an identity matrix. @@ -33,8 +38,12 @@ def test_reference_beam_mat(make_beam_rmat_impl, module_name): rmat = make_beam_rmat_impl(cnst.beam_vec, cnst.eta_vec) - assert_allclose(rmat, cnst.identity_3x3) + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + +# ------------------------------------------------------------------------------ + +# Test error conditions @all_impls def test_zero_beam_vec(make_beam_rmat_impl, module_name): @@ -51,6 +60,10 @@ def test_colinear_beam_eta_vec(make_beam_rmat_impl, module_name): make_beam_rmat_impl(cnst.beam_vec, cnst.beam_vec) +# ------------------------------------------------------------------------------ + +# Test orthonormal results + @all_impls def test_orthonormal_1(make_beam_rmat_impl, module_name): beam_vec = np.array([1.0, 2.0, 3.0]) @@ -66,7 +79,7 @@ def test_orthonormal_1(make_beam_rmat_impl, module_name): # dot(A, A.T) == Identity seems a good orthonormality check # Note: atol needed as rtol is not useful for '0.' entries. - assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=1e+10) + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=ATOL_IDENTITY) @all_impls @@ -82,7 +95,7 @@ def test_orthonormal_2(make_beam_rmat_impl, module_name): # dot(A, A.T) == Identity seems a good orthonormality check # Note: atol needed as rtol is not useful for '0.' entries. - assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=1e+10) + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=ATOL_IDENTITY) @all_impls @@ -95,4 +108,45 @@ def test_orthonormal_3(make_beam_rmat_impl, module_name): # dot(A, A.T) == Identity seems a good orthonormality check # Note: atol needed as rtol is not useful for '0.' entries. - assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=1e+10) + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, atol=ATOL_IDENTITY) + + +# ------------------------------------------------------------------------------ + +# Test strided inputs + +@all_impls +def test_strided_beam(make_beam_rmat_impl, module_name): + buff = np.zeros((3,2), order="C") + buff[:,:] = 42.0 # fill with some trash value + buff[:,0] = cnst.beam_vec # but set a strided vector to the valid value + + rmat = make_beam_rmat_impl(buff[:,0], cnst.eta_vec) + + # This should result in identity (see test_reference_beam_rmat) + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + + +@all_impls +def test_strided_eta(make_beam_rmat_impl, module_name): + buff = np.zeros((3,2), order="C") + buff[:,:] = 42.0 # fill with some trash value + buff[:,0] = cnst.eta_vec # but set a strided vector to the valid value + + rmat = make_beam_rmat_impl(cnst.beam_vec, buff[:,0]) + + # This should result in identity (see test_reference_beam_rmat) + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + + +@all_impls +def test_strided_beam_eta(make_beam_rmat_impl, module_name): + buff = np.zeros((3,4), order="C") + buff[:,:] = 42.0 + buff[:,0] = cnst.beam_vec + buff[:,2] = cnst.eta_vec + + rmat = make_beam_rmat_impl(buff[:,0], buff[:,2]) + # This should result in identity (see test_reference_beam_rmat) + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + From 74c5aa6e4aed3d3b5299cde39be5323e7cdb98bd Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 29 Oct 2018 18:29:20 +0100 Subject: [PATCH 23/30] first tests for make_rmat_of_expmap --- .../tests/test_make_rmat_of_expmap.py | 36 ++++++++++++++++--- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/hexrd/transforms/tests/test_make_rmat_of_expmap.py b/hexrd/transforms/tests/test_make_rmat_of_expmap.py index b33bacaa..dd36c955 100644 --- a/hexrd/transforms/tests/test_make_rmat_of_expmap.py +++ b/hexrd/transforms/tests/test_make_rmat_of_expmap.py @@ -7,8 +7,15 @@ from ..xf_capi import make_rmat_of_expmap as capi_make_rmat_of_expmap from ..xf_numba import make_rmat_of_expmap as numba_make_rmat_of_expmap +from ... import constants as cnst + +import numpy as np +from numpy.testing import assert_allclose + import pytest +ATOL_IDENTITY = 1e-10 + all_impls = pytest.mark.parametrize('make_rmat_of_expmap_impl, module_name', [(numpy_make_rmat_of_expmap, 'numpy'), (capi_make_rmat_of_expmap, 'capi'), @@ -17,10 +24,31 @@ ) + +# ------------------------------------------------------------------------------ + +# Test trivial case + @all_impls -def test_sample1(make_rmat_of_expmap_impl, module_name): - pass +def test_zero_expmap(make_rmat_of_expmap_impl, module_name): + exp_map = np.zeros((3,)) + + rmat = make_rmat_of_expmap_impl(exp_map) + + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + @all_impls -def test_sample2(make_rmat_of_expmap_impl, module_name): - pass +def test_2pi_expmap(make_rmat_of_expmap_impl, module_name): + """all this should result in identity - barring numerical error. + Note this goes via a different codepath as phi in the code is not 0.""" + + rmat = make_rmat_of_expmap_impl(np.array([2*np.pi, 0., 0.])) + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + + rmat = make_rmat_of_expmap_impl(np.array([0., 2*np.pi, 0.])) + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + + rmat = make_rmat_of_expmap_impl(np.array([0., 0.,2*np.pi])) + assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + From 94575e04cca9c84f06313bac8ff2a166ada7c949 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Mon, 29 Oct 2018 21:49:10 +0100 Subject: [PATCH 24/30] extra tests for make_rmat_of_expmap (orthonormal results). --- .../tests/test_make_rmat_of_expmap.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/hexrd/transforms/tests/test_make_rmat_of_expmap.py b/hexrd/transforms/tests/test_make_rmat_of_expmap.py index dd36c955..6fac54f7 100644 --- a/hexrd/transforms/tests/test_make_rmat_of_expmap.py +++ b/hexrd/transforms/tests/test_make_rmat_of_expmap.py @@ -52,3 +52,20 @@ def test_2pi_expmap(make_rmat_of_expmap_impl, module_name): rmat = make_rmat_of_expmap_impl(np.array([0., 0.,2*np.pi])) assert_allclose(rmat, cnst.identity_3x3, atol=ATOL_IDENTITY) + +# ------------------------------------------------------------------------------ + +# check that for some random inputs the resulting matrix is orthogonal + +@all_impls +def test_orthonormal(make_rmat_of_expmap_impl, module_name): + rmat = make_rmat_of_expmap_impl(np.array([42.0, 3., 32.5])) + # dot(A, A.T) == IDENTITY is a good orthonormality check + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, + atol=ATOL_IDENTITY) + + rmat = make_rmat_of_expmap_impl(np.array([-32.0, 0.0, 17.6])) + assert_allclose(np.dot(rmat, rmat.T), cnst.identity_3x3, + atol=ATOL_IDENTITY) + + From eb39871b4a95b2d673c86757384cc43b2caabb1d Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Tue, 30 Oct 2018 09:45:07 +0100 Subject: [PATCH 25/30] added strided test for make_rmat_of_expmap --- .../tests/test_make_rmat_of_expmap.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/hexrd/transforms/tests/test_make_rmat_of_expmap.py b/hexrd/transforms/tests/test_make_rmat_of_expmap.py index 6fac54f7..bc127be6 100644 --- a/hexrd/transforms/tests/test_make_rmat_of_expmap.py +++ b/hexrd/transforms/tests/test_make_rmat_of_expmap.py @@ -69,3 +69,20 @@ def test_orthonormal(make_rmat_of_expmap_impl, module_name): atol=ATOL_IDENTITY) + +# ------------------------------------------------------------------------------ + +# Test strided input +@all_impls +def test_strided(make_rmat_of_expmap_impl, module_name): + exp_map = np.array([42.0, 3., 32.5]) # A random expmap + + buff = np.zeros((3, 3), order='C') + buff[:,0] = exp_map[:] # assign the expmap to a column, so it is strided + + result_contiguous = make_rmat_of_expmap_impl(exp_map) + result_strided = make_rmat_of_expmap_impl(buff[:,0]) + + # in fact, a stricter equality check should work as well, + # but anyways... + assert_allclose(result_contiguous, result_strided) From 958e90634eb293a7f4986231238c07bf323c5b40 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Tue, 30 Oct 2018 10:41:46 +0100 Subject: [PATCH 26/30] added trivial test for unit_vector caught a bug in array version in numba. --- hexrd/transforms/tests/test_unit_vector.py | 23 +++++++++++++++++----- hexrd/transforms/xf_numba.py | 5 +++-- 2 files changed, 21 insertions(+), 7 deletions(-) diff --git a/hexrd/transforms/tests/test_unit_vector.py b/hexrd/transforms/tests/test_unit_vector.py index 47d85617..bc127228 100644 --- a/hexrd/transforms/tests/test_unit_vector.py +++ b/hexrd/transforms/tests/test_unit_vector.py @@ -7,6 +7,9 @@ from ..xf_capi import unit_vector as capi_unit_vector from ..xf_numba import unit_vector as numba_unit_vector +import numpy as np +from numpy.testing import assert_allclose + import pytest all_impls = pytest.mark.parametrize('unit_vector_impl, module_name', @@ -17,10 +20,20 @@ ) -@all_impls -def test_sample1(unit_vector_impl, module_name): - pass + +# ------------------------------------------------------------------------------ @all_impls -def test_sample2(unit_vector_impl, module_name): - pass +def test_trivial(unit_vector_impl, module_name): + # all vectors in eye(3) are already unit vectors + iden = np.eye(3) + + # check a vector at a time + assert_allclose(unit_vector_impl(iden[0]), iden[0]) + assert_allclose(unit_vector_impl(iden[1]), iden[1]) + assert_allclose(unit_vector_impl(iden[2]), iden[2]) + + # use the array version + assert_allclose(unit_vector_impl(iden), iden) + + diff --git a/hexrd/transforms/xf_numba.py b/hexrd/transforms/xf_numba.py index 9ec743bb..d722a86a 100644 --- a/hexrd/transforms/xf_numba.py +++ b/hexrd/transforms/xf_numba.py @@ -332,10 +332,10 @@ def _unit_vector_multi(a, out=None): nrm = np.sqrt(nrm) # prevent divide by zero if nrm > cnst.epsf: - for i in range(n): + for j in range(n): out[i, j] = a[i, j] / nrm else: - for i in range(n): + for j in range(n): out[i, j] = a[i, j] return out @@ -354,6 +354,7 @@ def row_norm(a): a = np.atleast_2d(a) return _row_norm(a) + @xf_api def unit_vector(vec_in): """ From 860e019d08b1939b832ffc26352a5e7e85691978 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Tue, 30 Oct 2018 10:52:58 +0100 Subject: [PATCH 27/30] added zero norm vector test in unit_vector --- hexrd/transforms/tests/test_unit_vector.py | 21 +++++++++++++++++++++ hexrd/transforms/transforms_definitions.py | 4 +++- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/hexrd/transforms/tests/test_unit_vector.py b/hexrd/transforms/tests/test_unit_vector.py index bc127228..02243f58 100644 --- a/hexrd/transforms/tests/test_unit_vector.py +++ b/hexrd/transforms/tests/test_unit_vector.py @@ -23,6 +23,8 @@ # ------------------------------------------------------------------------------ +# trivial tests + @all_impls def test_trivial(unit_vector_impl, module_name): # all vectors in eye(3) are already unit vectors @@ -37,3 +39,22 @@ def test_trivial(unit_vector_impl, module_name): assert_allclose(unit_vector_impl(iden), iden) +@all_impls +def test_zero(unit_vector_impl, module_name): + # When a zero vector is given, a potential "division by zero" happens. + # in this library, instead of trying to normalize a zero-norm vector + # (which would trigger the division by zero), the original vector is + # returned. + + # check vector + zero_vec = np.zeros((3,)) + assert_allclose(unit_vector_impl(zero_vec), zero_vec) + + # check array + zero_arr = np.zeros((3,3)) + assert_allclose(unit_vector_impl(zero_arr), zero_arr) + + # check mixed array + mix_arr = np.eye(3) + mix_arr[1,:] = 0.0 + assert_allclose(unit_vector_impl(mix_arr), mix_arr) diff --git a/hexrd/transforms/transforms_definitions.py b/hexrd/transforms/transforms_definitions.py index 30300bd5..a502f46b 100644 --- a/hexrd/transforms/transforms_definitions.py +++ b/hexrd/transforms/transforms_definitions.py @@ -363,7 +363,9 @@ def _signature(a): class DEF_unit_vector(DEF_Func): """ - normalize an array of row vectors (vstacked, axis=0) + Normalize an array of row vectors (vstacked, axis=0) + For vectors with (very close to zero) norm, the original + vector is returned. guaranteed to work for 1d and 2d arrays. """ From 6827224fa2a43328520867a08bdd011b41545a6f Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Tue, 30 Oct 2018 11:09:04 +0100 Subject: [PATCH 28/30] added unit vector test with some random vectors. --- hexrd/transforms/tests/test_unit_vector.py | 25 ++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/hexrd/transforms/tests/test_unit_vector.py b/hexrd/transforms/tests/test_unit_vector.py index 02243f58..6718926c 100644 --- a/hexrd/transforms/tests/test_unit_vector.py +++ b/hexrd/transforms/tests/test_unit_vector.py @@ -20,6 +20,16 @@ ) +def _get_random_vectors_array(): + # return a (n,3) array with some vectors and a (n) array with the expected + # result norms. + arr = np.array([[42.0, 0.0, 0.0, 1.0], + [12.0, 12.0, 12.0, 1.0], + [ 0.0, 0.0, 0.0, 0.0], + [ 0.7, -0.7, 0.0, 1.0], + [-0.0, -0.0, -0.0, 0.0]]) + return arr[:,0:3], arr[:,3] + # ------------------------------------------------------------------------------ @@ -58,3 +68,18 @@ def test_zero(unit_vector_impl, module_name): mix_arr = np.eye(3) mix_arr[1,:] = 0.0 assert_allclose(unit_vector_impl(mix_arr), mix_arr) + + +@all_impls +def test_random_vectors(unit_vector_impl, module_name): + # test for some random vectors. The test just checks that the norm of the + # the resulting vector is as expected. + + vecs, expected_norm = _get_random_vectors_array() + + # element by element + for i in range(len(vecs)): + assert_allclose(np.linalg.norm(unit_vector_impl(vecs[i])), expected_norm[i]) + + # all in a row + assert_allclose(np.linalg.norm(unit_vector_impl(vecs), axis=1), expected_norm) From 1a8ad9ac672d4dc8cabd6960219bb7081cf81337 Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Tue, 30 Oct 2018 11:21:23 +0100 Subject: [PATCH 29/30] added unit vector test with strided input vectors. --- hexrd/transforms/tests/test_unit_vector.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/hexrd/transforms/tests/test_unit_vector.py b/hexrd/transforms/tests/test_unit_vector.py index 6718926c..68ba7136 100644 --- a/hexrd/transforms/tests/test_unit_vector.py +++ b/hexrd/transforms/tests/test_unit_vector.py @@ -83,3 +83,21 @@ def test_random_vectors(unit_vector_impl, module_name): # all in a row assert_allclose(np.linalg.norm(unit_vector_impl(vecs), axis=1), expected_norm) + + +# ------------------------------------------------------------------------------ + +# check input ordering + +@all_impls +def test_strided_inputs(unit_vector_impl, module_name): + vecs, expected_norm = _get_random_vectors_array() + + vecs_f = np.asfortranarray(vecs) + + # element by element + for i in range(len(vecs_f)): + assert_allclose(np.linalg.norm(unit_vector_impl(vecs_f[i])), expected_norm[i]) + + # all in a row + assert_allclose(np.linalg.norm(unit_vector_impl(vecs_f), axis=1), expected_norm) From aeab000627b238e92c5122579cce5b354827c48e Mon Sep 17 00:00:00 2001 From: Oscar Villellas Date: Tue, 30 Oct 2018 18:08:19 +0100 Subject: [PATCH 30/30] Added a README for transform tests So that it serves as documentation for the submodule testing as well as guidelines for test writers. --- hexrd/transforms/tests/README.md | 91 ++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 hexrd/transforms/tests/README.md diff --git a/hexrd/transforms/tests/README.md b/hexrd/transforms/tests/README.md new file mode 100644 index 00000000..dbc5bd8e --- /dev/null +++ b/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 ` 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_.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): + +``` + +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.