From 32009f305d6e28905785e499e17d58eb38663333 Mon Sep 17 00:00:00 2001 From: Kevin Lewis <55998034+kevindlewis23@users.noreply.github.com> Date: Mon, 12 Aug 2024 11:16:49 -0400 Subject: [PATCH 01/16] Crystallography Refactored (#701) * Reformat utils * Remove unused _convert_angles function * Remove check for numba install * Clean up unused variables, add typing * Remove unused _coo_build_window function * Add type hints to function defs * Remove unnecessary variables / passes * Provide deprecation warnings to the user. * Tuple -> tuple * List -> list, Dict -> dict * Change deprecated to accept str instead of func * Remove unnecessary "save" func from etaOmeMaps * Add deprecation date * Format with Black * Clean up name mangling, bad returns * Move to f-strings * Clean up properties * Delete some unused functions * Remove makeNew function, which is not a useful method * Some docstrings, committing to swtich branches * Some docstrings written, properties now using annotations * Typing and cleanup * Remove unused functions/args * Added docstrings * Update crystallography.py * Push before pull * Crystallography done * PEP8 * Union instead of | * Copy plane data --------- Co-authored-by: Zack Singer --- hexrd/findorientations.py | 6 +- hexrd/instrument/detector.py | 9 +- hexrd/material/crystallography.py | 1702 ++++++++++++++--------------- tests/planedata/test_exclusion.py | 11 +- tests/planedata/test_misc.py | 3 + tests/planedata/test_with_data.py | 3 + 6 files changed, 846 insertions(+), 888 deletions(-) diff --git a/hexrd/findorientations.py b/hexrd/findorientations.py index b04fe4fd3..aa24aeb80 100755 --- a/hexrd/findorientations.py +++ b/hexrd/findorientations.py @@ -88,7 +88,7 @@ def generate_orientation_fibers(cfg, eta_ome): pd = eta_ome.planeData tTh = pd.getTTh() bMat = pd.latVecOps['B'] - csym = pd.getLaueGroup() + csym = pd.laueGroup # !!! changed recently where iHKLList are now master hklIDs pd_hkl_ids = eta_ome.iHKLList[seed_hkl_ids] @@ -885,8 +885,8 @@ def find_orientations(cfg, logger.info("\tmean reflections per grain: %d", mean_rpg) logger.info("\tneighborhood size: %d", min_samples) - qbar, _ = run_cluster( - completeness, qfib, plane_data.getQSym(), cfg, + qbar, cl = run_cluster( + completeness, qfib, plane_data.q_sym, cfg, min_samples=min_samples, compl_thresh=compl_thresh, radius=cl_radius) diff --git a/hexrd/instrument/detector.py b/hexrd/instrument/detector.py index 9071a4896..227f20370 100644 --- a/hexrd/instrument/detector.py +++ b/hexrd/instrument/detector.py @@ -1131,14 +1131,7 @@ def make_powder_rings( # in case you want to give it tth angles directly if isinstance(pd, PlaneData): - # Okay, we have a PlaneData object - try: - pd = pd.makeNew() # make a copy to munge - except TypeError: - # !!! have some other object here, - # likely a dummy plane data object of sorts - raise - + pd = PlaneData(None, pd) if delta_tth is not None: pd.tThWidth = np.radians(delta_tth) else: diff --git a/hexrd/material/crystallography.py b/hexrd/material/crystallography.py index 76d74fd77..13148d295 100644 --- a/hexrd/material/crystallography.py +++ b/hexrd/material/crystallography.py @@ -28,16 +28,21 @@ import re import copy from math import pi +from typing import Optional, Union, Dict, List, Tuple import numpy as np -import csv -import os +from hexrd.material.unitcell import unitcell +from hexrd.deprecation import deprecated from hexrd import constants from hexrd.matrixutil import unitVector -from hexrd.rotations import \ - rotMatOfExpMap, mapAngle, applySym, \ - ltypeOfLaueGroup, quatOfLaueGroup +from hexrd.rotations import ( + rotMatOfExpMap, + mapAngle, + applySym, + ltypeOfLaueGroup, + quatOfLaueGroup, +) from hexrd.transforms import xfcapi from hexrd import valunits from hexrd.valunits import toFloat @@ -51,24 +56,60 @@ outputDegrees_bak = outputDegrees -def hklToStr(x): - return re.sub(r'[\[\]\(\)\{\},]', '', str(x)) +def hklToStr(hkl: np.ndarray) -> str: + """ + Converts hkl representation to a string. + + Parameters + ---------- + hkl : np.ndarray + 3 element list of h, k, and l values (Miller indices). + + Returns + ------- + str + Space-separated string representation of h, k, and l values. + + """ + return re.sub(r'[\[\]\(\)\{\},]', '', str(hkl)) + + +def tempSetOutputDegrees(val: bool) -> None: + """ + Set the global outputDegrees flag temporarily. Can be reverted with + revertOutputDegrees(). + + Parameters + ---------- + val : bool + True to output angles in degrees, False to output angles in radians. + Returns + ------- + None -def tempSetOutputDegrees(val): + """ global outputDegrees, outputDegrees_bak outputDegrees_bak = outputDegrees outputDegrees = val - return -def revertOutputDegrees(): +def revertOutputDegrees() -> None: + """ + Revert the effect of tempSetOutputDegrees(), resetting the outputDegrees + flag to its previous value (True to output in degrees, False for radians). + + Returns + ------- + None + """ global outputDegrees, outputDegrees_bak outputDegrees = outputDegrees_bak - return -def cosineXform(a, b, c): +def cosineXform( + a: np.ndarray, b: np.ndarray, c: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: """ Spherical trig transform to take alpha, beta, gamma to expressions for cos(alpha*). See ref below. @@ -77,67 +118,59 @@ def cosineXform(a, b, c): the relations between direct and reciprocal lattice quantities''. Acta Cryst. (1968), A24, 247--248 + Parameters + ---------- + a : np.ndarray + List of alpha angle values (radians). + b : np.ndarray + List of beta angle values (radians). + c : np.ndarray + List of gamma angle values (radians). + + Returns + ------- + np.ndarray + List of cos(alpha*) values. + np.ndarray + List of sin(alpha*) values. + """ - cosar = (np.cos(b)*np.cos(c) - np.cos(a)) / (np.sin(b)*np.sin(c)) + cosar = (np.cos(b) * np.cos(c) - np.cos(a)) / (np.sin(b) * np.sin(c)) sinar = np.sqrt(1 - cosar**2) return cosar, sinar -def processWavelength(arg): +def processWavelength(arg: Union[valunits.valWUnit, float]) -> float: """ Convert an energy value to a wavelength. If argument has units of length or energy, will convert to globally specified unit type for wavelength (dUnit). If argument is a scalar, assumed input units are keV. """ - if hasattr(arg, 'getVal'): + if isinstance(arg, valunits.valWUnit): + # arg is a valunits.valWUnit object if arg.isLength(): - retval = arg.getVal(dUnit) + return arg.getVal(dUnit) elif arg.isEnergy(): e = arg.getVal('keV') - retval = valunits.valWUnit( + return valunits.valWUnit( 'wavelength', 'length', constants.keVToAngstrom(e), 'angstrom' ).getVal(dUnit) else: - raise RuntimeError('do not know what to do with '+str(arg)) + raise RuntimeError('do not know what to do with ' + str(arg)) else: # !!! assuming arg is in keV - retval = valunits.valWUnit( + return valunits.valWUnit( 'wavelength', 'length', constants.keVToAngstrom(arg), 'angstrom' ).getVal(dUnit) - return retval - - -def latticeParameters(lvec): - """ - Generates direct and reciprocal lattice vector components in a - crystal-relative RHON basis, X. The convention for fixing X to the - lattice is such that a || x1 and c* || x3, where a and c* are - direct and reciprocal lattice vectors, respectively. - """ - lnorm = np.sqrt(np.sum(lvec**2, 0)) - - a = lnorm[0] - b = lnorm[1] - c = lnorm[2] - - ahat = lvec[:, 0]/a - bhat = lvec[:, 1]/b - chat = lvec[:, 2]/c - gama = np.arccos(np.dot(ahat, bhat)) - beta = np.arccos(np.dot(ahat, chat)) - alfa = np.arccos(np.dot(bhat, chat)) - if outputDegrees: - gama = r2d*gama - beta = r2d*beta - alfa = r2d*alfa - - return [a, b, c, alfa, beta, gama] - - -def latticePlanes(hkls, lparms, - ltype='cubic', wavelength=1.54059292, strainMag=None): +def latticePlanes( + hkls: np.ndarray, + lparms: np.ndarray, + ltype: Optional[str] = 'cubic', + wavelength: Optional[float] = 1.54059292, + strainMag: Optional[float] = None, +) -> Dict[str, np.ndarray]: """ Generates lattice plane data in the direct lattice for a given set of Miller indices. Vector components are written in the @@ -158,9 +191,9 @@ def latticePlanes(hkls, lparms, 2) lparms (1 x m float list) is the array of lattice parameters, where m depends on the symmetry group (see below). - 3) The following optional keyword arguments are recognized: + The following optional arguments are recognized: - *) ltype=(string) is a string representing the symmetry type of + 3) ltype=(string) is a string representing the symmetry type of the implied Laue group. The 11 available choices are shown below. The default value is 'cubic'. Note that each group expects a lattice parameter array of the indicated length @@ -177,11 +210,11 @@ def latticePlanes(hkls, lparms, 'monoclinic' a, b, c, beta (in degrees) 'triclinic' a, b, c, alpha, beta, gamma (in degrees) - *) wavelength= is a value represented the wavelength in + 4) wavelength= is a value represented the wavelength in Angstroms to calculate bragg angles for. The default value is for Cu K-alpha radiation (1.54059292 Angstrom) - *) strainMag=None + 5) strainMag=None OUTPUTS: @@ -194,7 +227,7 @@ def latticePlanes(hkls, lparms, dspacings (n, ) double array array of the d-spacings for each {hkl} - 2thetas (n, ) double array array of the Bragg angles for + tThetas (n, ) double array array of the Bragg angles for each {hkl} relative to the specified wavelength @@ -216,8 +249,9 @@ def latticePlanes(hkls, lparms, """ location = 'latticePlanes' - assert hkls.shape[0] == 3, \ - "hkls aren't column vectors in call to '%s'!" % location + assert ( + hkls.shape[0] == 3 + ), f"hkls aren't column vectors in call to '{location}'!" tag = ltype wlen = wavelength @@ -231,43 +265,46 @@ def latticePlanes(hkls, lparms, # magnitudes d = 1 / np.sqrt(np.sum(G**2, 0)) - aconv = 1. + aconv = 1.0 if outputDegrees: aconv = r2d # two thetas - sth = wlen / 2. / d - mask = (np.abs(sth) < 1.) + sth = wlen / 2.0 / d + mask = np.abs(sth) < 1.0 tth = np.zeros(sth.shape) tth[~mask] = np.nan - tth[mask] = aconv * 2. * np.arcsin(sth[mask]) + tth[mask] = aconv * 2.0 * np.arcsin(sth[mask]) - p = dict(normals=unitVector(G), - dspacings=d, - tThetas=tth) + p = dict(normals=unitVector(G), dspacings=d, tThetas=tth) if strainMag is not None: p['tThetasLo'] = np.zeros(sth.shape) p['tThetasHi'] = np.zeros(sth.shape) - mask = ( - (np.abs(wlen / 2. / (d * (1. + strainMag))) < 1.) & - (np.abs(wlen / 2. / (d * (1. - strainMag))) < 1.) + mask = (np.abs(wlen / 2.0 / (d * (1.0 + strainMag))) < 1.0) & ( + np.abs(wlen / 2.0 / (d * (1.0 - strainMag))) < 1.0 ) p['tThetasLo'][~mask] = np.nan p['tThetasHi'][~mask] = np.nan - p['tThetasLo'][mask] = \ - aconv * 2 * np.arcsin(wlen/2./(d[mask]*(1. + strainMag))) - p['tThetasHi'][mask] = \ - aconv * 2 * np.arcsin(wlen/2./(d[mask]*(1. - strainMag))) + p['tThetasLo'][mask] = ( + aconv * 2 * np.arcsin(wlen / 2.0 / (d[mask] * (1.0 + strainMag))) + ) + p['tThetasHi'][mask] = ( + aconv * 2 * np.arcsin(wlen / 2.0 / (d[mask] * (1.0 - strainMag))) + ) return p -def latticeVectors(lparms, tag='cubic', radians=False, debug=False): +def latticeVectors( + lparms: np.ndarray, + tag: Optional[str] = 'cubic', + radians: Optional[bool] = False, +) -> Dict[str, Union[np.ndarray, float]]: """ Generates direct and reciprocal lattice vector components in a crystal-relative RHON basis, X. The convention for fixing X to the @@ -283,7 +320,7 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): 1) lparms (1 x n float list) is the array of lattice parameters, where n depends on the symmetry group (see below). - 2) symTag (string) is a case-insensitive string representing the + 2) tag (string) is a case-insensitive string representing the symmetry type of the implied Laue group. The 11 available choices are shown below. The default value is 'cubic'. Note that each group expects a lattice parameter array of the indicated length @@ -300,6 +337,11 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): 'monoclinic' a, b, c, beta (in degrees) 'triclinic' a, b, c, alpha, beta, gamma (in degrees) + The following optional arguments are recognized: + + 3) radians= is a boolean flag indicating usage of radians rather + than degrees, defaults to false. + OUTPUTS: 1) lattice is a dictionary containing the following keys/items: @@ -380,60 +422,64 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): 'tetragonal', 'orthorhombic', 'monoclinic', - 'triclinic' - ] + 'triclinic', + ] if radians: - aconv = 1. + aconv = 1.0 else: - aconv = pi/180. # degToRad - deg90 = pi/2. - deg120 = 2.*pi/3. + aconv = pi / 180.0 # degToRad + deg90 = pi / 2.0 + deg120 = 2.0 * pi / 3.0 # if tag == lattStrings[0]: # cubic - cellparms = np.r_[np.tile(lparms[0], (3,)), deg90*np.ones((3,))] + cellparms = np.r_[np.tile(lparms[0], (3,)), deg90 * np.ones((3,))] elif tag == lattStrings[1] or tag == lattStrings[2]: # hexagonal | trigonal (hex indices) - cellparms = np.r_[lparms[0], lparms[0], lparms[1], - deg90, deg90, deg120] + cellparms = np.r_[ + lparms[0], lparms[0], lparms[1], deg90, deg90, deg120 + ] elif tag == lattStrings[3]: # rhombohedral - cellparms = np.r_[np.tile(lparms[0], (3,)), - np.tile(aconv*lparms[1], (3,))] + cellparms = np.r_[ + np.tile(lparms[0], (3,)), np.tile(aconv * lparms[1], (3,)) + ] elif tag == lattStrings[4]: # tetragonal - cellparms = np.r_[lparms[0], lparms[0], lparms[1], - deg90, deg90, deg90] + cellparms = np.r_[lparms[0], lparms[0], lparms[1], deg90, deg90, deg90] elif tag == lattStrings[5]: # orthorhombic - cellparms = np.r_[lparms[0], lparms[1], lparms[2], - deg90, deg90, deg90] + cellparms = np.r_[lparms[0], lparms[1], lparms[2], deg90, deg90, deg90] elif tag == lattStrings[6]: # monoclinic - cellparms = np.r_[lparms[0], lparms[1], lparms[2], - deg90, aconv*lparms[3], deg90] + cellparms = np.r_[ + lparms[0], lparms[1], lparms[2], deg90, aconv * lparms[3], deg90 + ] elif tag == lattStrings[7]: # triclinic - # FIXME: fixed DP 2/24/16 - cellparms = np.r_[lparms[0], lparms[1], lparms[2], - aconv*lparms[3], aconv*lparms[4], aconv*lparms[5]] + cellparms = np.r_[ + lparms[0], + lparms[1], + lparms[2], + aconv * lparms[3], + aconv * lparms[4], + aconv * lparms[5], + ] else: - raise RuntimeError('lattice tag \'%s\' is not recognized' % (tag)) - - if debug: - print((str(cellparms[0:3]) + ' ' + str(r2d*cellparms[3:6]))) - alfa = cellparms[3] - beta = cellparms[4] - gama = cellparms[5] + raise RuntimeError(f'lattice tag "{tag}" is not recognized') - cosalfar, sinalfar = cosineXform(alfa, beta, gama) + alpha, beta, gamma = cellparms[3:6] + cosalfar, sinalfar = cosineXform(alpha, beta, gamma) - a = cellparms[0]*np.r_[1, 0, 0] - b = cellparms[1]*np.r_[np.cos(gama), np.sin(gama), 0] - c = cellparms[2]*np.r_[np.cos(beta), - -cosalfar*np.sin(beta), - sinalfar*np.sin(beta)] + a = cellparms[0] * np.r_[1, 0, 0] + b = cellparms[1] * np.r_[np.cos(gamma), np.sin(gamma), 0] + c = ( + cellparms[2] + * np.r_[ + np.cos(beta), -cosalfar * np.sin(beta), sinalfar * np.sin(beta) + ] + ) ad = np.sqrt(np.sum(a**2)) bd = np.sqrt(np.sum(b**2)) @@ -446,200 +492,62 @@ def latticeVectors(lparms, tag='cubic', radians=False, debug=False): F = np.c_[a, b, c] # Reciprocal lattice vectors - astar = np.cross(b, c)/V - bstar = np.cross(c, a)/V - cstar = np.cross(a, b)/V + astar = np.cross(b, c) / V + bstar = np.cross(c, a) / V + cstar = np.cross(a, b) / V # and parameters ar = np.sqrt(np.sum(astar**2)) br = np.sqrt(np.sum(bstar**2)) cr = np.sqrt(np.sum(cstar**2)) - alfar = np.arccos(np.dot(bstar, cstar)/br/cr) - betar = np.arccos(np.dot(cstar, astar)/cr/ar) - gamar = np.arccos(np.dot(astar, bstar)/ar/br) + alfar = np.arccos(np.dot(bstar, cstar) / br / cr) + betar = np.arccos(np.dot(cstar, astar) / cr / ar) + gamar = np.arccos(np.dot(astar, bstar) / ar / br) # B takes components in the reciprocal lattice to X B = np.c_[astar, bstar, cstar] cosalfar2, sinalfar2 = cosineXform(alfar, betar, gamar) - afable = ar*np.r_[1, 0, 0] - bfable = br*np.r_[np.cos(gamar), np.sin(gamar), 0] - cfable = cr*np.r_[np.cos(betar), - -cosalfar2*np.sin(betar), - sinalfar2*np.sin(betar)] + afable = ar * np.r_[1, 0, 0] + bfable = br * np.r_[np.cos(gamar), np.sin(gamar), 0] + cfable = ( + cr + * np.r_[ + np.cos(betar), + -cosalfar2 * np.sin(betar), + sinalfar2 * np.sin(betar), + ] + ) BR = np.c_[afable, bfable, cfable] U0 = np.dot(B, np.linalg.inv(BR)) if outputDegrees: - dparms = np.r_[ad, bd, cd, r2d*np.r_[alfa, beta, gama]] - rparms = np.r_[ar, br, cr, r2d*np.r_[alfar, betar, gamar]] + dparms = np.r_[ad, bd, cd, r2d * np.r_[alpha, beta, gamma]] + rparms = np.r_[ar, br, cr, r2d * np.r_[alfar, betar, gamar]] else: - dparms = np.r_[ad, bd, cd, np.r_[alfa, beta, gama]] + dparms = np.r_[ad, bd, cd, np.r_[alpha, beta, gamma]] rparms = np.r_[ar, br, cr, np.r_[alfar, betar, gamar]] - L = {'F': F, - 'B': B, - 'BR': BR, - 'U0': U0, - 'vol': V, - 'dparms': dparms, - 'rparms': rparms} - - return L - - -def hexagonalIndicesFromRhombohedral(hkl): - """ - converts rhombohedral hkl to hexagonal indices - """ - HKL = np.zeros((3, hkl.shape[1]), dtype='int') - - HKL[0, :] = hkl[0, :] - hkl[1, :] - HKL[1, :] = hkl[1, :] - hkl[2, :] - HKL[2, :] = hkl[0, :] + hkl[1, :] + hkl[2, :] - - return HKL - - -def rhombohedralIndicesFromHexagonal(HKL): - """ - converts hexagonal hkl to rhombohedral indices - """ - hkl = np.zeros((3, HKL.shape[1]), dtype='int') - - hkl[0, :] = 2 * HKL[0, :] + HKL[1, :] + HKL[2, :] - hkl[1, :] = -HKL[0, :] + HKL[1, :] + HKL[2, :] - hkl[2, :] = -HKL[0, :] - 2 * HKL[1, :] + HKL[2, :] - - hkl = hkl / 3. - return hkl - - -def rhombohedralParametersFromHexagonal(a_h, c_h): - """ - converts hexagonal lattice parameters (a, c) to rhombohedral - lattice parameters (a, alpha) - """ - a_r = np.sqrt(3 * a_h**2 + c_h**2) / 3. - alfa_r = 2 * np.arcsin(3. / (2 * np.sqrt(3 + (c_h / a_h)**2))) - if outputDegrees: - alfa_r = r2d * alfa_r - return a_r, alfa_r - - -def convert_Miller_direction_to_cartesian(uvw, a=1., c=1., normalize=False): - """ - Converts 3-index hexagonal Miller direction indices to components in the - crystal reference frame. - - Parameters - ---------- - uvw : array_like - The (n, 3) array of 3-index hexagonal indices to convert. - a : scalar, optional - The `a` lattice parameter. The default value is 1. - c : scalar, optional - The `c` lattice parameter. The default value is 1. - normalize : bool, optional - Flag for whether or not to normalize output vectors - - Returns - ------- - numpy.ndarray - The (n, 3) array of cartesian components associated with the input - direction indices. - - Notes - ----- - 1) The [uv.w] the Miller-Bravais convention is in the hexagonal basis - {a1, a2, a3, c}. The basis for the output, {o1, o2, o3}, is - chosen such that - - o1 || a1 - o3 || c - o2 = o3 ^ o1 - - """ - u, v, w = np.atleast_2d(uvw).T - retval = np.vstack([1.5*u*a, sqrt3by2*a*(2*v + u), w*c]) - if normalize: - return unitVector(retval).T - else: - return retval.T - - -def convert_Miller_direction_to_MillerBravias(uvw, suppress_redundant=True): - """ - Converts 3-index hexagonal Miller direction indices to 4-index - Miller-Bravais direction indices. - - Parameters - ---------- - uvw : array_like - The (n, 3) array of 3-index hexagonal Miller indices to convert. - suppress_redundant : bool, optional - Flag to suppress the redundant 3rd index. The default is True. - - Returns - ------- - numpy.ndarray - The (n, 3) or (n, 4) array -- depending on kwarg -- of Miller-Bravis - components associated with the input Miller direction indices. - - Notes - ----- - * NOT for plane normals!!! - - """ - u, v, w = np.atleast_2d(uvw).T - retval = np.vstack([(2*u - v)/3, (2*v - u)/3, w]).T - rem = np.vstack([np.mod(np.tile(i[0], 2), i[1:]) for i in retval]) - rem[abs(rem) < epsf] = np.nan - lcm = np.nanmin(rem, axis=1) - lcm[np.isnan(lcm)] = 1 - retval = retval / np.tile(lcm, (3, 1)).T - if suppress_redundant: - return retval - else: - t = np.atleast_2d(1 - np.sum(retval[:2], axis=1)).T - return np.hstack([retval[:, :2], t, np.atleast_2d(retval[:, 2]).T]) - - -def convert_MillerBravias_direction_to_Miller(UVW): - """ - Converts 4-index hexagonal Miller-Bravais direction indices to - 3-index Miller direction indices. - - Parameters - ---------- - UVW : array_like - The (n, 3) array of **non-redundant** Miller-Bravais direction indices - to convert. - - Returns - ------- - numpy.ndarray - The (n, 3) array of Miller direction indices associated with the - input Miller-Bravais indices. - - Notes - ----- - * NOT for plane normals!!! - - """ - U, V, W = np.atleast_2d(UVW).T - return np.vstack([2*U + V, 2*V + U, W]) + return { + 'F': F, + 'B': B, + 'BR': BR, + 'U0': U0, + 'vol': V, + 'dparms': dparms, + 'rparms': rparms, + } class PlaneData(object): """ Careful with ordering: Outputs are ordered by the 2-theta for the - hkl unless you get self.__hkls directly, and this order can change + hkl unless you get self._hkls directly, and this order can change with changes in lattice parameters (lparms); setting and getting exclusions works on the current hkl ordering, not the original - ordering (in self.__hkls), but exclusions are stored in the + ordering (in self._hkls), but exclusions are stored in the original ordering in case the hkl ordering does change with lattice parameters @@ -648,77 +556,94 @@ class PlaneData(object): tThWidth """ - def __init__(self, - hkls, - *args, - **kwargs): + def __init__(self, hkls: Optional[np.ndarray], *args, **kwargs) -> None: + """ + Constructor for PlaneData + + Parameters + ---------- + hkls : np.ndarray + Miller indices to be used in the plane data. Can be None if + args is another PlaneData object + + *args + Unnamed arguments. Could be in the format of `lparms, laueGroup, + wavelength, strainMag`, or just a `PlaneData` object. + + **kwargs + Valid keyword arguments include: + - doTThSort + - exclusions + - tThMax + - tThWidth + """ + self._doTThSort = True + self._exclusions = None + self._tThMax = None - self.phaseID = None - self.__doTThSort = True - self.__exclusions = None - self.__tThMax = None - # if len(args) == 4: lparms, laueGroup, wavelength, strainMag = args tThWidth = None - self.__wavelength = processWavelength(wavelength) - self.__lparms = self.__parseLParms(lparms) - elif len(args) == 1 and hasattr(args[0], 'getParams'): + self._wavelength = processWavelength(wavelength) + self._lparms = self._parseLParms(lparms) + elif len(args) == 1 and isinstance(args[0], PlaneData): other = args[0] - lparms, laueGroup, wavelength, strainMag, tThWidth = \ + lparms, laueGroup, wavelength, strainMag, tThWidth = ( other.getParams() - self.__wavelength = wavelength - self.__lparms = lparms - self.phaseID = other.phaseID - self.__doTThSort = other.__doTThSort - self.__exclusions = other.__exclusions - self.__tThMax = other.__tThMax + ) + self._wavelength = wavelength + self._lparms = lparms + self._doTThSort = other._doTThSort + self._exclusions = other._exclusions + self._tThMax = other._tThMax if hkls is None: - hkls = other.__hkls + hkls = other._hkls else: - raise NotImplementedError('args : '+str(args)) + raise NotImplementedError(f'args : {args}') - self.__laueGroup = laueGroup - self.__qsym = quatOfLaueGroup(self.__laueGroup) - self.__hkls = copy.deepcopy(hkls) - self.__strainMag = strainMag - self.__structFact = np.ones(self.__hkls.shape[1]) + self._laueGroup = laueGroup + self._hkls = copy.deepcopy(hkls) + self._strainMag = strainMag + self._structFact = np.ones(self._hkls.shape[1]) self.tThWidth = tThWidth # ... need to implement tThMin too - if 'phaseID' in kwargs: - self.phaseID = kwargs.pop('phaseID') if 'doTThSort' in kwargs: - self.__doTThSort = kwargs.pop('doTThSort') + self._doTThSort = kwargs.pop('doTThSort') if 'exclusions' in kwargs: - self.__exclusions = kwargs.pop('exclusions') + self._exclusions = kwargs.pop('exclusions') if 'tThMax' in kwargs: - self.__tThMax = toFloat(kwargs.pop('tThMax'), 'radians') + self._tThMax = toFloat(kwargs.pop('tThMax'), 'radians') if 'tThWidth' in kwargs: self.tThWidth = kwargs.pop('tThWidth') if len(kwargs) > 0: - raise RuntimeError('have unparsed keyword arguments with keys: ' - + str(list(kwargs.keys()))) + raise RuntimeError( + f'have unparsed keyword arguments with keys: {kwargs.keys()}' + ) # This is only used to calculate the structure factor if invalidated - self.__unitcell = None - - self.__calc() - - return - - def __calc(self): - symmGroup = ltypeOfLaueGroup(self.__laueGroup) - latPlaneData, latVecOps, hklDataList = PlaneData.makePlaneData( - self.__hkls, self.__lparms, self.__qsym, - symmGroup, self.__strainMag, self.wavelength) + self._unitcell: unitcell = None + + self._calc() + + def _calc(self): + symmGroup = ltypeOfLaueGroup(self._laueGroup) + self._q_sym = quatOfLaueGroup(self._laueGroup) + _, latVecOps, hklDataList = PlaneData.makePlaneData( + self._hkls, + self._lparms, + self._q_sym, + symmGroup, + self._strainMag, + self.wavelength, + ) 'sort by tTheta' tThs = np.array( [hklDataList[iHKL]['tTheta'] for iHKL in range(len(hklDataList))] ) - if self.__doTThSort: - # sorted hkl -> __hkl - # __hkl -> sorted hkl + if self._doTThSort: + # sorted hkl -> _hkl + # _hkl -> sorted hkl self.tThSort = np.argsort(tThs) self.tThSortInv = np.empty(len(hklDataList), dtype=int) self.tThSortInv[self.tThSort] = np.arange(len(hklDataList)) @@ -729,101 +654,131 @@ def __calc(self): self.hklDataList = hklDataList self._latVecOps = latVecOps self.nHKLs = len(self.getHKLs()) - return def __str__(self): s = '========== plane data ==========\n' s += 'lattice parameters:\n ' + str(self.lparms) + '\n' - s += 'two theta width: (%s)\n' % str(self.tThWidth) - s += 'strain magnitude: (%s)\n' % str(self.strainMag) - s += 'beam energy (%s)\n' % str(self.wavelength) + s += f'two theta width: ({str(self.tThWidth)})\n' + s += f'strain magnitude: ({str(self.strainMag)})\n' + s += f'beam energy ({str(self.wavelength)})\n' s += 'hkls: (%d)\n' % self.nHKLs s += str(self.getHKLs()) return s - def getNHKLs(self): - return self.nHKLs + def getParams(self): + """ + Getter for the parameters of the plane data. + + Returns + ------- + tuple + The parameters of the plane data. In the order of + _lparams, _laueGroup, _wavelength, _strainMag, tThWidth - def getPhaseID(self): - 'may return None if not set' - return self.phaseID + """ + return ( + self._lparms, + self._laueGroup, + self._wavelength, + self._strainMag, + self.tThWidth, + ) - def getParams(self): - return (self.__lparms, self.__laueGroup, self.__wavelength, - self.__strainMag, self.tThWidth) + def getNhklRef(self) -> int: + """ + Get the total number of hkl's in the plane data, not ignoring + ones that are excluded in exclusions. - def getNhklRef(self): - 'does not use exclusions or the like' - retval = len(self.hklDataList) - return retval + Returns + ------- + int + The total number of hkl's in the plane data. + """ + return len(self.hklDataList) - def get_hkls(self): + @property + def hkls(self) -> np.ndarray: """ - do not do return self.__hkls, as everywhere else hkls are returned - in 2-theta order; transpose is to comply with lparm convention + hStacked Hkls of the plane data (Miller indices). """ return self.getHKLs().T - def set_hkls(self, hkls): - raise RuntimeError('for now, not allowing hkls to be reset') - # self.__exclusions = None - # self.__hkls = hkls - # self.__calc() - return - hkls = property(get_hkls, set_hkls, None) + @hkls.setter + def hkls(self, hkls): + raise NotImplementedError('for now, not allowing hkls to be reset') + + @property + def tThMax(self) -> Optional[float]: + """ + Maximum 2-theta value of the plane data. - def get_tThMax(self): - return self.__tThMax + float or None + """ + return self._tThMax - def set_tThMax(self, tThMax): - self.__tThMax = toFloat(tThMax, 'radians') - # self.__calc() # no need to redo calc for tThMax - return - tThMax = property(get_tThMax, set_tThMax, None) + @tThMax.setter + def tThMax(self, t_th_max: Union[float, valunits.valWUnit]) -> None: + self._tThMax = toFloat(t_th_max, 'radians') - def get_exclusions(self): + @property + def exclusions(self) -> np.ndarray: + """ + Excluded HKL's the plane data. + + Set as type np.ndarray, as a mask of length getNhklRef(), a list of + indices to be excluded, or a list of ranges of indices. + + Read as a mask of length getNhklRef(). + """ retval = np.zeros(self.getNhklRef(), dtype=bool) - if self.__exclusions is not None: + if self._exclusions is not None: # report in current hkl ordering - retval[:] = self.__exclusions[self.tThSortInv] - if self.__tThMax is not None: + retval[:] = self._exclusions[self.tThSortInv] + if self._tThMax is not None: for iHKLr, hklData in enumerate(self.hklDataList): - if hklData['tTheta'] > self.__tThMax: + if hklData['tTheta'] > self._tThMax: retval[iHKLr] = True return retval - def set_exclusions(self, exclusions): + @exclusions.setter + def exclusions(self, new_exclusions: Optional[np.ndarray]) -> None: excl = np.zeros(len(self.hklDataList), dtype=bool) - if exclusions is not None: - exclusions = np.atleast_1d(exclusions) + if new_exclusions is not None: + exclusions = np.atleast_1d(new_exclusions) if len(exclusions) == len(self.hklDataList): - assert exclusions.dtype == 'bool', \ - 'exclusions should be bool if full length' - # convert from current hkl ordering to __hkl ordering + assert ( + exclusions.dtype == 'bool' + ), 'Exclusions should be bool if full length' + # convert from current hkl ordering to _hkl ordering excl[:] = exclusions[self.tThSort] else: if len(exclusions.shape) == 1: # treat exclusions as indices excl[self.tThSort[exclusions]] = True elif len(exclusions.shape) == 2: - raise NotImplementedError( - 'have not yet coded treating exclusions as ranges' - ) + # treat exclusions as ranges of indices + for r in exclusions: + excl[self.tThSort[r[0]:r[1]]] = True else: raise RuntimeError( - 'do not now what to do with exclusions with shape ' - + str(exclusions.shape) + f'Unclear behavior for shape {exclusions.shape}' ) - self.__exclusions = excl - self.nHKLs = np.sum(np.logical_not(self.__exclusions)) - return - exclusions = property(get_exclusions, set_exclusions, None) + self._exclusions = excl + self.nHKLs = np.sum(np.logical_not(self._exclusions)) def exclude( - self, dmin=None, dmax=None, tthmin=None, tthmax=None, - sfacmin=None, sfacmax=None, pintmin=None, pintmax=None - ): - """Set exclusions according to various parameters + self, + dmin: Optional[float] = None, + dmax: Optional[float] = None, + tthmin: Optional[float] = None, + tthmax: Optional[float] = None, + sfacmin: Optional[float] = None, + sfacmax: Optional[float] = None, + pintmin: Optional[float] = None, + pintmax: Optional[float] = None, + ) -> None: + """ + Set exclusions according to various parameters Any hkl with a value below any min or above any max will be excluded. So to be included, an hkl needs to have values between the min and max @@ -856,43 +811,42 @@ def exclude( if (dmin is not None) or (dmax is not None): d = np.array(self.getPlaneSpacings()) - if (dmin is not None): + if dmin is not None: excl[d < dmin] = True - if (dmax is not None): + if dmax is not None: excl[d > dmax] = True if (tthmin is not None) or (tthmax is not None): tth = self.getTTh() - if (tthmin is not None): + if tthmin is not None: excl[tth < tthmin] = True - if (tthmax is not None): + if tthmax is not None: excl[tth > tthmax] = True if (sfacmin is not None) or (sfacmax is not None): sfac = self.structFact - sfac = sfac/sfac.max() - if (sfacmin is not None): + sfac = sfac / sfac.max() + if sfacmin is not None: excl[sfac < sfacmin] = True - if (sfacmax is not None): + if sfacmax is not None: excl[sfac > sfacmax] = True if (pintmin is not None) or (pintmax is not None): pint = self.powder_intensity - pint = pint/pint.max() - if (pintmin is not None): + pint = pint / pint.max() + if pintmin is not None: excl[pint < pintmin] = True - if (pintmax is not None): + if pintmax is not None: excl[pint > pintmax] = True self.exclusions = excl - def get_lparms(self): - return self.__lparms - - def __parseLParms(self, lparms): + def _parseLParms( + self, lparms: List[Union[valunits.valWUnit, float]] + ) -> List[float]: lparmsDUnit = [] for lparmThis in lparms: - if hasattr(lparmThis, 'getVal'): + if isinstance(lparmThis, valunits.valWUnit): if lparmThis.isLength(): lparmsDUnit.append(lparmThis.getVal(dUnit)) elif lparmThis.isAngle(): @@ -901,70 +855,104 @@ def __parseLParms(self, lparms): lparmsDUnit.append(lparmThis.getVal('degrees')) else: raise RuntimeError( - 'do not know what to do with ' - + str(lparmThis) + f'Do not know what to do with {lparmThis}' ) else: lparmsDUnit.append(lparmThis) return lparmsDUnit - def set_lparms(self, lparms): - self.__lparms = self.__parseLParms(lparms) - self.__calc() - return - lparms = property(get_lparms, set_lparms, None) + @property + def lparms(self) -> List[float]: + """ + Lattice parameters of the plane data. + + Can be set as a List[float | valWUnit], but will be converted to + List[float]. + """ + return self._lparms - def get_strainMag(self): - return self.__strainMag + @lparms.setter + def lparms(self, lparms: List[Union[valunits.valWUnit, float]]) -> None: + self._lparms = self._parseLParms(lparms) + self._calc() + + @property + def strainMag(self) -> Optional[float]: + """ + Strain magnitude of the plane data. + + float or None + """ + return self._strainMag - def set_strainMag(self, strainMag): - self.__strainMag = strainMag + @strainMag.setter + def strainMag(self, strain_mag: float) -> None: + self._strainMag = strain_mag self.tThWidth = None - self.__calc() - return - strainMag = property(get_strainMag, set_strainMag, None) + self._calc() - def get_wavelength(self): - return self.__wavelength + @property + def wavelength(self) -> float: + """ + Wavelength of the plane data. + + Set as float or valWUnit. - def set_wavelength(self, wavelength): + Read as float + """ + return self._wavelength + + @wavelength.setter + def wavelength(self, wavelength: Union[float, valunits.valWUnit]) -> None: wavelength = processWavelength(wavelength) - if np.isclose(self.__wavelength, wavelength): - # Do not re-compute if it is almost the same + # Do not re-compute if it is almost the same + if np.isclose(self._wavelength, wavelength): return - self.__wavelength = wavelength - self.__calc() - - wavelength = property(get_wavelength, set_wavelength, None) + self._wavelength = wavelength + self._calc() - def invalidate_structure_factor(self, unitcell): - # It can be expensive to compute the structure factor, so provide the - # option to just invalidate it, while providing a unit cell, so that - # it can be lazily computed from the unit cell. - self.__structFact = None + def invalidate_structure_factor(self, ucell: unitcell) -> None: + """ + It can be expensive to compute the structure factor + This method just invalidates it, providing a unit cell, + so that it can be lazily computed from the unit cell. + + Parameters: + ----------- + unitcell : unitcell + The unit cell to be used to compute the structure factor + """ + self._structFact = None self._hedm_intensity = None self._powder_intensity = None - self.__unitcell = unitcell + self._unitcell = ucell def _compute_sf_if_needed(self): any_invalid = ( - self.__structFact is None or - self._hedm_intensity is None or - self._powder_intensity is None + self._structFact is None + or self._hedm_intensity is None + or self._powder_intensity is None ) - if any_invalid and self.__unitcell is not None: + if any_invalid and self._unitcell is not None: # Compute the structure factor first. # This can be expensive to do, so we lazily compute it when needed. hkls = self.getHKLs(allHKLs=True) - self.set_structFact(self.__unitcell.CalcXRSF(hkls)) + self.structFact = self._unitcell.CalcXRSF(hkls) + + @property + def structFact(self) -> np.ndarray: + """ + Structure factors for each hkl. - def get_structFact(self): + np.ndarray + """ self._compute_sf_if_needed() - return self.__structFact[~self.exclusions] + return self._structFact[~self.exclusions] - def set_structFact(self, structFact): - self.__structFact = structFact + @structFact.setter + def structFact(self, structFact: np.ndarray) -> None: + self._structFact = structFact multiplicity = self.getMultiplicity(allHKLs=True) tth = self.getTTh(allHKLs=True) @@ -981,33 +969,70 @@ def set_structFact(self, structFact): self._hedm_intensity = hedm_intensity self._powder_intensity = powderI - structFact = property(get_structFact, set_structFact, None) - @property - def powder_intensity(self): + def powder_intensity(self) -> np.ndarray: + """ + Powder intensity for each hkl. + """ self._compute_sf_if_needed() return self._powder_intensity[~self.exclusions] @property - def hedm_intensity(self): + def hedm_intensity(self) -> np.ndarray: + """ + HEDM (high energy x-ray diffraction microscopy) intensity for each hkl. + """ self._compute_sf_if_needed() return self._hedm_intensity[~self.exclusions] @staticmethod - def makePlaneData(hkls, lparms, qsym, symmGroup, strainMag, wavelength): + def makePlaneData( + hkls: np.ndarray, + lparms: np.ndarray, + qsym: np.ndarray, + symmGroup, + strainMag, + wavelength, + ) -> Tuple[ + Dict[str, np.ndarray], Dict[str, Union[np.ndarray, float]], List[Dict] + ]: """ - hkls : need to work with crystallography.latticePlanes - lparms : need to work with crystallography.latticePlanes - laueGroup : see symmetry module - wavelength : wavelength - strainMag : swag of strian magnitudes + Generate lattice plane data from inputs. + + Parameters: + ----------- + hkls: np.ndarray + Miller indices, as in crystallography.latticePlanes + lparms: np.ndarray + Lattice parameters, as in crystallography.latticePlanes + qsym: np.ndarray + (4, n) containing quaternions of symmetry + symmGroup: str + Tag for the symmetry (Laue) group of the lattice. Can generate from + ltypeOfLaueGroup + strainMag: float + Swag of strain magnitudes + wavelength: float + Wavelength + + Returns: + ------- + dict: + Dictionary containing lattice plane data + dict: + Dictionary containing lattice vector operators + list: + List of dictionaries, each containing the data for one hkl """ tempSetOutputDegrees(False) - latPlaneData = latticePlanes(hkls, lparms, - ltype=symmGroup, - strainMag=strainMag, - wavelength=wavelength) + latPlaneData = latticePlanes( + hkls, + lparms, + ltype=symmGroup, + strainMag=strainMag, + wavelength=wavelength, + ) latVecOps = latticeVectors(lparms, symmGroup) @@ -1027,113 +1052,128 @@ def makePlaneData(hkls, lparms, qsym, symmGroup, strainMag, wavelength): np.dot(latVecOps['B'], hkls[:, iHKL].reshape(3, 1)), qsym, csFlag=True, - cullPM=False) + cullPM=False, + ) # check for +/- in symmetry group latPlnNrmlsM = applySym( np.dot(latVecOps['B'], hkls[:, iHKL].reshape(3, 1)), qsym, csFlag=False, - cullPM=False) + cullPM=False, + ) csRefl = latPlnNrmls.shape[1] == latPlnNrmlsM.shape[1] # added this so that I retain the actual symmetric # integer hkls as well symHKLs = np.array( - np.round( - np.dot(latVecOps['F'].T, latPlnNrmls) - ), dtype='int' + np.round(np.dot(latVecOps['F'].T, latPlnNrmls)), dtype='int' ) hklDataList.append( - dict(hklID=iHKL, - hkl=hkls[:, iHKL], - tTheta=latPlaneData['tThetas'][iHKL], - dSpacings=latPlaneData['dspacings'][iHKL], - tThetaLo=latPlaneData['tThetasLo'][iHKL], - tThetaHi=latPlaneData['tThetasHi'][iHKL], - latPlnNrmls=unitVector(latPlnNrmls), - symHKLs=symHKLs, - centrosym=csRefl - ) + dict( + hklID=iHKL, + hkl=hkls[:, iHKL], + tTheta=latPlaneData['tThetas'][iHKL], + dSpacings=latPlaneData['dspacings'][iHKL], + tThetaLo=latPlaneData['tThetasLo'][iHKL], + tThetaHi=latPlaneData['tThetasHi'][iHKL], + latPlnNrmls=unitVector(latPlnNrmls), + symHKLs=symHKLs, + centrosym=csRefl, ) + ) revertOutputDegrees() return latPlaneData, latVecOps, hklDataList - def getLatticeType(self): - """This is the lattice type""" - return ltypeOfLaueGroup(self.__laueGroup) - - def getLaueGroup(self): - """This is the Schoenflies tag""" - return self.__laueGroup + @property + def laueGroup(self) -> str: + """ + This is the Schoenflies tag, describing symmetry group of the lattice. + Note that setting this with incompatible lattice parameters will + cause an error. If changing both, use set_laue_and_lparms. - def setLaueGroup(self, laueGroup): - self.__laueGroup = laueGroup - self.__calc() + str + """ + return self._laueGroup - laueGroup = property(getLaueGroup, setLaueGroup, None) + @laueGroup.setter + def laueGroup(self, laueGroup: str) -> None: + self._laueGroup = laueGroup + self._calc() - def set_laue_and_lparms(self, laueGroup, lparms): - """Set the Laue group and lattice parameters simultaneously + def set_laue_and_lparms( + self, laueGroup: str, lparms: List[Union[valunits.valWUnit, float]] + ) -> None: + """ + Set the Laue group and lattice parameters simultaneously When the Laue group changes, the lattice parameters may be - incompatible, and cause an error in self.__calc(). This function + incompatible, and cause an error in self._calc(). This function allows us to update both the Laue group and lattice parameters simultaneously to avoid this issue. + + Parameters: + ----------- + laueGroup : str + The symmetry (Laue) group to be set + lparms : List[valunits.valWUnit | float] + Lattice parameters to be set """ - self.__laueGroup = laueGroup - self.__lparms = self.__parseLParms(lparms) - self.__calc() + self._laueGroup = laueGroup + self._lparms = self._parseLParms(lparms) + self._calc() - def getQSym(self): - return self.__qsym # rotations.quatOfLaueGroup(self.__laueGroup) + @property + def q_sym(self) -> np.ndarray: + """ + Quaternions of symmetry for each hkl, generated from the Laue group + + np.ndarray((4, n)) + """ + return self._q_sym # rotations.quatOfLaueGroup(self._laueGroup) - def getPlaneSpacings(self): + def getPlaneSpacings(self) -> List[float]: """ - gets plane spacings + Plane spacings for each hkl. + + Returns: + ------- + List[float] + List of plane spacings for each hkl """ dspacings = [] for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue dspacings.append(hklData['dSpacings']) return dspacings - def getPlaneNormals(self): - """ - gets both +(hkl) and -(hkl) normals - """ - plnNrmls = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - plnNrmls.append(hklData['latPlnNrmls']) - return plnNrmls - @property - def latVecOps(self): + def latVecOps(self) -> Dict[str, Union[np.ndarray, float]]: """ gets lattice vector operators as a new (deepcopy) + + Returns: + ------- + Dict[str, np.ndarray | float] + Dictionary containing lattice vector operators """ return copy.deepcopy(self._latVecOps) - def __thisHKL(self, iHKLr): - retval = True + def _thisHKL(self, iHKLr: int) -> bool: hklData = self.hklDataList[iHKLr] - if self.__exclusions is not None: - if self.__exclusions[self.tThSortInv[iHKLr]]: - retval = False - if self.__tThMax is not None: - # FIXME: check for nans here??? - if hklData['tTheta'] > self.__tThMax \ - or np.isnan(hklData['tTheta']): - retval = False - return retval - - def __getTThRange(self, iHKLr): + if self._exclusions is not None: + if self._exclusions[self.tThSortInv[iHKLr]]: + return False + if self._tThMax is not None: + if hklData['tTheta'] > self._tThMax or np.isnan(hklData['tTheta']): + return False + return True + + def _getTThRange(self, iHKLr: int) -> Tuple[float, float]: hklData = self.hklDataList[iHKLr] if self.tThWidth is not None: # tThHi-tThLo < self.tThWidth tTh = hklData['tTheta'] @@ -1144,40 +1184,58 @@ def __getTThRange(self, iHKLr): tThLo = hklData['tThetaLo'] return (tThLo, tThHi) - def getTThRanges(self, strainMag=None, lparms=None): + def getTThRanges(self, strainMag: Optional[float] = None) -> np.ndarray: """ - Return 2-theta ranges for included hkls + Get the 2-theta ranges for included hkls + + Parameters: + ----------- + strainMag : Optional[float] + Optional swag of strain magnitude - return array is n x 2 + Returns: + ------- + np.ndarray: + hstacked array of hstacked tThLo and tThHi for each hkl (n x 2) """ - if lparms is None: - tThRanges = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - # tThRanges.append([hklData['tThetaLo'], hklData['tThetaHi']]) - if strainMag is None: - tThRanges.append(self.__getTThRange(iHKLr)) - else: - hklData = self.hklDataList[iHKLr] - d = hklData['dSpacings'] - tThLo = 2.0 * np.arcsin( - self.__wavelength / 2. / (d*(1.+strainMag)) - ) - tThHi = 2.0 * np.arcsin( - self.__wavelength / 2. / (d*(1.-strainMag)) - ) - tThRanges.append((tThLo, tThHi)) - else: - new = self.__class__(self.__hkls, self) - new.lparms = lparms - tThRanges = new.getTThRanges(strainMag=strainMag) + tThRanges = [] + for iHKLr, hklData in enumerate(self.hklDataList): + if not self._thisHKL(iHKLr): + continue + if strainMag is None: + tThRanges.append(self._getTThRange(iHKLr)) + else: + hklData = self.hklDataList[iHKLr] + d = hklData['dSpacings'] + tThLo = 2.0 * np.arcsin( + self._wavelength / 2.0 / (d * (1.0 + strainMag)) + ) + tThHi = 2.0 * np.arcsin( + self._wavelength / 2.0 / (d * (1.0 - strainMag)) + ) + tThRanges.append((tThLo, tThHi)) return np.array(tThRanges) - def getMergedRanges(self, cullDupl=False): + def getMergedRanges( + self, cullDupl: Optional[bool] = False + ) -> Tuple[List[List[int]], List[List[float]]]: """ - return indices and ranges for specified planeData, merging where + Return indices and ranges for specified planeData, merging where there is overlap based on the tThWidth and line positions + + Parameters: + ----------- + cullDupl : (optional) bool + If True, cull duplicate 2-theta values (within sqrt_epsf). Defaults + to False. + + Returns: + -------- + List[List[int]] + List of indices for each merged range + + List[List[float]] + List of merged ranges, (n x 2) """ tThs = self.getTTh() tThRanges = self.getTThRanges() @@ -1191,11 +1249,11 @@ def getMergedRanges(self, cullDupl=False): mergedRanges = [] hklsCur = [] tThLoIdx = 0 - tThHiCur = 0. + tThHiCur = 0.0 for iHKL, nonoverlapNext in enumerate(nonoverlapNexts): tThHi = tThRanges[iHKL, -1] if not nonoverlapNext: - if cullDupl and abs(tThs[iHKL] - tThs[iHKL+1]) < sqrt_epsf: + if cullDupl and abs(tThs[iHKL] - tThs[iHKL + 1]) < sqrt_epsf: continue else: hklsCur.append(iHKL) @@ -1209,71 +1267,55 @@ def getMergedRanges(self, cullDupl=False): hklsCur = [] return iHKLLists, mergedRanges - def makeNew(self): - new = self.__class__(None, self) - return new - - def getTTh(self, lparms=None, allHKLs=False): - if lparms is None: - tTh = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not allHKLs: - if not self.__thisHKL(iHKLr): - continue - tTh.append(hklData['tTheta']) - else: - tTh.append(hklData['tTheta']) - else: - new = self.makeNew() - new.lparms = lparms - tTh = new.getTTh() - return np.array(tTh) - - def getDD_tThs_lparms(self): - """ - derivatives of tThs with respect to lattice parameters; - have not yet done coding for analytic derivatives, just wimp out - and finite difference + def getTTh(self, allHKLs: Optional[bool] = False) -> np.ndarray: """ - pert = 1.0e-5 # assume they are all around unity - pertInv = 1.0/pert + Get the 2-theta values for each hkl. - lparmsRef = copy.deepcopy(self.__lparms) - tThRef = self.getTTh() - ddtTh = np.empty((len(tThRef), len(lparmsRef))) + Parameters: + ----------- + allHKLs : (optional) bool + If True, return all 2-theta values, even if they are excluded in + the current planeData. Default is False. - for iLparm in range(len(lparmsRef)): - self.__lparms = copy.deepcopy(lparmsRef) - self.__lparms[iLparm] += pert - self.__calc() - - iTTh = 0 - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - ddtTh[iTTh, iLparm] = \ - np.r_[hklData['tTheta'] - tThRef[iTTh]]*pertInv - iTTh += 1 + Returns: + ------- + np.ndarray + Array of 2-theta values for each hkl + """ + tTh = [] + for iHKLr, hklData in enumerate(self.hklDataList): + if not allHKLs and not self._thisHKL(iHKLr): + continue + tTh.append(hklData['tTheta']) + return np.array(tTh) - 'restore' - self.__lparms = lparmsRef - self.__calc() + def getMultiplicity(self, allHKLs: Optional[bool] = False) -> np.ndarray: + """ + Get the multiplicity for each hkl (number of symHKLs). - return ddtTh + Paramters: + ---------- + allHKLs : (optional) bool + If True, return all multiplicities, even if they are excluded in + the current planeData. Defaults to false. - def getMultiplicity(self, allHKLs=False): + Returns + ------- + np.ndarray + Array of multiplicities for each hkl + """ # ... JVB: is this incorrect? multip = [] for iHKLr, hklData in enumerate(self.hklDataList): - if not allHKLs: - if not self.__thisHKL(iHKLr): - continue - multip.append(hklData['symHKLs'].shape[1]) - else: + if allHKLs or self._thisHKL(iHKLr): multip.append(hklData['symHKLs'].shape[1]) return np.array(multip) - def getHKLID(self, hkl, master=False): + def getHKLID( + self, + hkl: Union[int, Tuple[int, int, int], np.ndarray], + master: Optional[bool] = False, + ) -> Union[List[int], int]: """ Return the unique ID of a list of hkls. @@ -1293,7 +1335,7 @@ def getHKLID(self, hkl, master=False): Returns ------- - retval : list + hkl_ids : list The list of requested hklID values associate with the input. Notes @@ -1305,23 +1347,24 @@ def getHKLID(self, hkl, master=False): 2020-05-21 (JVB) -- modified to handle all symmetric equavlent reprs. """ if hasattr(hkl, '__setitem__'): # tuple does not have __setitem__ - if hasattr(hkl, 'shape'): + if isinstance(hkl, np.ndarray): # if is ndarray, assume is 3xN - # retval = list(map(self.__getHKLID, hkl.T)) - retval = [self.__getHKLID(x, master=master) for x in hkl.T] + return [self._getHKLID(x, master=master) for x in hkl.T] else: - # retval = list(map(self.__getHKLID, hkl)) - retval = [self.__getHKLID(x, master=master) for x in hkl] + return [self._getHKLID(x, master=master) for x in hkl] else: - retval = self.__getHKLID(hkl, master=master) - return retval + return self._getHKLID(hkl, master=master) - def __getHKLID(self, hkl, master=False): + def _getHKLID( + self, + hkl: Union[int, Tuple[int, int, int], np.ndarray], + master: Optional[bool] = False, + ) -> int: """ for hkl that is a tuple, return externally visible hkl index """ if isinstance(hkl, int): - retval = hkl + return hkl else: hklList = self.getSymHKLs() # !!! list, reduced by exclusions intl_hklIDs = np.asarray([i['hklID'] for i in self.hklDataList]) @@ -1332,14 +1375,13 @@ def __getHKLID(self, hkl, master=False): for thisHKL in symHKLs.T: dHKLInv[tuple(thisHKL)] = idx try: - retval = dHKLInv[tuple(hkl)] - except(KeyError): + return dHKLInv[tuple(hkl)] + except KeyError: raise RuntimeError( f"hkl '{tuple(hkl)}' is not present in this material!" ) - return retval - def getHKLs(self, *hkl_ids, **kwargs): + def getHKLs(self, *hkl_ids: int, **kwargs) -> Union[List[str], np.ndarray]: """ Returns the powder HKLs subject to specified options. @@ -1366,7 +1408,7 @@ def getHKLs(self, *hkl_ids, **kwargs): Returns ------- - retval : list | numpy.ndarray + hkls : list | numpy.ndarray Either a list of hkls as strings (if asStr=True) or a vstacked array of hkls. @@ -1396,10 +1438,10 @@ def getHKLs(self, *hkl_ids, **kwargs): if len(hkl_ids) == 0: for iHKLr, hklData in enumerate(self.hklDataList): if not opts['allHKLs']: - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue if opts['thisTTh'] is not None: - tThLo, tThHi = self.__getTThRange(iHKLr) + tThLo, tThHi = self._getTThRange(iHKLr) if opts['thisTTh'] < tThHi and opts['thisTTh'] > tThLo: hkls.append(hklData['hkl']) else: @@ -1415,7 +1457,7 @@ def getHKLs(self, *hkl_ids, **kwargs): # find ordinal index of current hklID try: idx[i] = int(np.where(all_hkl_ids == hkl_id)[0]) - except(TypeError): + except TypeError: raise RuntimeError( f"Requested hklID '{hkl_id}'is invalid!" ) @@ -1427,88 +1469,104 @@ def getHKLs(self, *hkl_ids, **kwargs): # handle output kwarg if opts['asStr']: - retval = list(map(hklToStr, np.array(hkls))) + return list(map(hklToStr, np.array(hkls))) else: - retval = np.array(hkls) - return retval - - def getSymHKLs(self, asStr=False, withID=False, indices=None): + return np.array(hkls) + + def getSymHKLs( + self, + asStr: Optional[bool] = False, + withID: Optional[bool] = False, + indices: Optional[List[int]] = None, + ) -> Union[List[List[str]], List[np.ndarray]]: """ - new function that returns all symmetric hkls + Return all symmetry HKLs. + + Parameters + ---------- + asStr : bool, optional + If True, return the symmetry HKLs as strings. The default is False. + withID : bool, optional + If True, return the symmetry HKLs with the hklID. The default is + False. Does nothing if asStr is True. + indices : list[inr], optional + Optional list of indices of hkls to include. + + Returns + ------- + sym_hkls : list list of strings, or list of numpy.ndarray + List of symmetry HKLs for each HKL, either as strings or as a + vstacked array. """ - retval = [] - iRetval = 0 + sym_hkls = [] + hkl_index = 0 if indices is not None: indB = np.zeros(self.nHKLs, dtype=bool) indB[np.array(indices)] = True else: indB = np.ones(self.nHKLs, dtype=bool) for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue - if indB[iRetval]: + if indB[hkl_index]: hkls = hklData['symHKLs'] if asStr: - retval.append(list(map(hklToStr, np.array(hkls).T))) + sym_hkls.append(list(map(hklToStr, np.array(hkls).T))) elif withID: - retval.append( + sym_hkls.append( np.vstack( - [np.tile(hklData['hklID'], (1, hkls.shape[1])), - hkls] + [ + np.tile(hklData['hklID'], (1, hkls.shape[1])), + hkls, + ] ) ) else: - retval.append(np.array(hkls)) - iRetval += 1 - return retval - - def getCentroSymHKLs(self): - retval = [] - for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): - continue - retval.append(hklData['centrosym']) - return retval - - def makeTheseScatteringVectors(self, hklList, rMat, - bMat=None, wavelength=None, chiTilt=None): - iHKLList = np.atleast_1d(self.getHKLID(hklList)) - fHKLs = np.hstack(self.getSymHKLs(indices=iHKLList)) - if bMat is None: - bMat = self._latVecOps['B'] - if wavelength is None: - wavelength = self.__wavelength - retval = PlaneData.makeScatteringVectors( - fHKLs, rMat, bMat, wavelength, chiTilt=chiTilt - ) - return retval - - def makeAllScatteringVectors(self, rMat, - bMat=None, wavelength=None, chiTilt=None): - fHKLs = np.hstack(self.getSymHKLs()) - if bMat is None: - bMat = self._latVecOps['B'] - if wavelength is None: - wavelength = self.__wavelength - retval = PlaneData.makeScatteringVectors( - fHKLs, rMat, bMat, wavelength, chiTilt=chiTilt - ) - return retval + sym_hkls.append(np.array(hkls)) + hkl_index += 1 + return sym_hkls @staticmethod - def makeScatteringVectors(hkls, rMat_c, bMat, wavelength, chiTilt=None): + def makeScatteringVectors( + hkls: np.ndarray, + rMat_c: np.ndarray, + bMat: np.ndarray, + wavelength: float, + chiTilt: Optional[float] = None, + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ Static method for calculating g-vectors and scattering vector angles for specified hkls, subject to the bragg conditions specified by lattice vectors, orientation matrix, and wavelength + Parameters + ---------- + hkls : np.ndarray + (3, n) array of hkls. + rMat_c : np.ndarray + (3, 3) rotation matrix from the crystal to the sample frame. + bMat : np.ndarray, optional + (3, 3) COB from reciprocal lattice frame to the crystal frame. + wavelength : float + xray wavelength in Angstroms. + chiTilt : float, optional + 0 <= chiTilt <= 90 degrees, defaults to 0 + + Returns + ------- + gVec_s : np.ndarray + (3, n) array of g-vectors (reciprocal lattice) in the sample frame. + oangs0 : np.ndarray + (3, n) array containing the feasible (2-theta, eta, ome) triplets + for each input hkl (first solution) + oangs1 : np.ndarray + (3, n) array containing the feasible (2-theta, eta, ome) triplets + for each input hkl (second solution) + FIXME: must do testing on strained bMat """ # arg munging - if chiTilt is None: - chi = 0. - else: - chi = float(chiTilt) + chi = float(chiTilt) if chiTilt is not None else 0.0 rMat_c = rMat_c.squeeze() # these are the reciprocal lattice vectors in the SAMPLE FRAME @@ -1517,9 +1575,9 @@ def makeScatteringVectors(hkls, rMat_c, bMat, wavelength, chiTilt=None): # strained [a, b, c] in the CRYSTAL FRAME gVec_s = np.dot(rMat_c, np.dot(bMat, hkls)) - dim0, nRefl = gVec_s.shape - assert dim0 == 3, \ - "Looks like something is wrong with your lattice plane normals" + dim0 = gVec_s.shape[0] + if dim0 != 3: + raise ValueError(f'Number of lattice plane normal dims is {dim0}') # call model from transforms now oangs0, oangs1 = xfcapi.oscill_angles_of_hkls( @@ -1528,7 +1586,12 @@ def makeScatteringVectors(hkls, rMat_c, bMat, wavelength, chiTilt=None): return gVec_s, oangs0.T, oangs1.T - def __makeScatteringVectors(self, rMat, bMat=None, chiTilt=None): + def _makeScatteringVectors( + self, + rMat: np.ndarray, + bMat: Optional[np.ndarray] = None, + chiTilt: Optional[float] = None, + ) -> Tuple[List[np.ndarray], List[np.ndarray], List[np.ndarray]]: """ modeled after QFromU.m """ @@ -1540,11 +1603,14 @@ def __makeScatteringVectors(self, rMat, bMat=None, chiTilt=None): Qs_ang0 = [] Qs_ang1 = [] for iHKLr, hklData in enumerate(self.hklDataList): - if not self.__thisHKL(iHKLr): + if not self._thisHKL(iHKLr): continue thisQs, thisAng0, thisAng1 = PlaneData.makeScatteringVectors( - hklData['symHKLs'], rMat, bMat, - self.__wavelength, chiTilt=chiTilt + hklData['symHKLs'], + rMat, + bMat, + self._wavelength, + chiTilt=chiTilt, ) Qs_vec.append(thisQs) Qs_ang0.append(thisAng0) @@ -1552,70 +1618,30 @@ def __makeScatteringVectors(self, rMat, bMat=None, chiTilt=None): return Qs_vec, Qs_ang0, Qs_ang1 - def calcStructFactor(self, atominfo): - """ - Calculates unit cell structure factors as a function of hkl - - USAGE: - - FSquared = calcStructFactor(atominfo,hkls,B) - - INPUTS: - - 1) atominfo (m x 1 float ndarray) the first threee columns of the - matrix contain fractional atom positions [uvw] of atoms in the unit - cell. The last column contains the number of electrons for a given atom - - 2) hkls (3 x n float ndarray) is the array of Miller indices for - the planes of interest. The vectors are assumed to be - concatenated along the 1-axis (horizontal) - - 3) B (3 x 3 float ndarray) is a matrix of reciprocal lattice basis - vectors,where each column contains a reciprocal lattice basis vector - ({g}=[B]*{hkl}) - + # OLD DEPRECATED PLANE_DATA STUFF ==================================== + @deprecated(new_func="len(self.hkls.T)", removal_date="2025-08-01") + def getNHKLs(self): + return len(self.getHKLs()) - OUTPUTS: + @deprecated(new_func="self.exclusions", removal_date="2025-08-01") + def get_exclusions(self): + return self.exclusions - 1) FSquared (n x 1 float ndarray) array of structure factors, - one for each hkl passed into the function - """ - r = atominfo[:, 0:3] - elecNum = atominfo[:, 3] - hkls = self.hkls - B = self.latVecOps['B'] - sinThOverLamdaList, ffDataList = LoadFormFactorData() - FSquared = np.zeros(hkls.shape[1]) - - for jj in np.arange(0, hkls.shape[1]): - # ???: probably have other functions for this - # Calculate G for each hkl - # Calculate magnitude of G for each hkl - G = hkls[0, jj]*B[:, 0] + hkls[1, jj]*B[:, 1] + hkls[2, jj]*B[:, 2] - magG = np.sqrt(G[0]**2 + G[1]**2 + G[2]**2) - - # Begin calculating form factor - F = 0 - for ii in np.arange(0, r.shape[0]): - ff = RetrieveAtomicFormFactor( - elecNum[ii], magG, sinThOverLamdaList, ffDataList - ) - exparg = complex( - 0., - 2.*np.pi*(hkls[0, jj]*r[ii, 0] - + hkls[1, jj]*r[ii, 1] - + hkls[2, jj]*r[ii, 2]) - ) - F = F + ff*np.exp(exparg) + @deprecated(new_func="self.exclusions=...", removal_date="2025-08-01") + def set_exclusions(self, exclusions): + self.exclusions = exclusions - """ - F = sum_atoms(ff(Q)*e^(2*pi*i(hu+kv+lw))) - """ - FSquared[jj] = np.real(F*np.conj(F)) + @deprecated(new_func="rotations.ltypeOfLaueGroup(self.laueGroup)", + removal_date="2025-08-01") + def getLatticeType(self): + return ltypeOfLaueGroup(self.laueGroup) - return FSquared + @deprecated(new_func="self.q_sym", removal_date="2025-08-01") + def getQSym(self): + return self.q_sym +@deprecated(removal_date='2025-01-01') def getFriedelPair(tth0, eta0, *ome0, **kwargs): """ Get the diffractometer angular coordinates in degrees for @@ -1692,22 +1718,12 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): dispFlag = False fableFlag = False chi = None - c1 = 1. - c2 = pi/180. - - # cast to arrays (in case they aren't) - if np.isscalar(eta0): - eta0 = [eta0] - - if np.isscalar(tth0): - tth0 = [tth0] - - if np.isscalar(ome0): - ome0 = [ome0] + c1 = 1.0 + c2 = pi / 180.0 - eta0 = np.asarray(eta0) - tth0 = np.asarray(tth0) - ome0 = np.asarray(ome0) + eta0 = np.atleast_1d(eta0) + tth0 = np.atleast_1d(tth0) + ome0 = np.atleast_1d(ome0) if eta0.ndim != 1: raise RuntimeError('azimuthal input must be 1-D') @@ -1719,25 +1735,24 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): else: if len(tth0) != npts: if len(tth0) == 1: - tth0 = tth0*np.ones(npts) + tth0 *= np.ones(npts) elif npts == 1: npts = len(tth0) - eta0 = eta0*np.ones(npts) + eta0 *= np.ones(npts) else: raise RuntimeError( 'the azimuthal and Bragg angle inputs are inconsistent' ) if len(ome0) == 0: - ome0 = np.zeros(npts) # dummy ome0 + ome0 = np.zeros(npts) # dummy ome0 elif len(ome0) == 1 and npts > 1: - ome0 = ome0*np.ones(npts) + ome0 *= np.ones(npts) else: if len(ome0) != npts: raise RuntimeError( 'your oscialltion angle input is inconsistent; ' - + 'it has length %d while it should be %d' - % (len(ome0), npts) + + f'it has length {len(ome0)} while it should be {npts}' ) # keyword args processing @@ -1752,8 +1767,8 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): fableFlag = True elif argkeys[i] == 'units': if kwargs[argkeys[i]] == 'radians': - c1 = 180./pi - c2 = 1. + c1 = 180.0 / pi + c2 = 1.0 elif argkeys[i] == 'chiTilt': if kwargs[argkeys[i]] is not None: chi = kwargs[argkeys[i]] @@ -1767,17 +1782,17 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): # mapped eta input # - in DEGREES, thanks to c1 - eta0 = mapAngle(c1*eta0, [-180, 180], units='degrees') + eta0 = mapAngle(c1 * eta0, [-180, 180], units='degrees') if fableFlag: eta0 = 90 - eta0 # must put args into RADIANS # - eta0 is in DEGREES, # - the others are in whatever was entered, hence c2 - eta0 = d2r*eta0 - tht0 = c2*tth0/2 + eta0 = d2r * eta0 + tht0 = c2 * tth0 / 2 if chi is not None: - chi = c2*chi + chi = c2 * chi else: chi = 0 @@ -1789,15 +1804,15 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): = sin(theta) - sin(chi)sin(eta)cos(theta) - Identity: a sin x + b cos x = sqrt(a**2 + b**2) sin (x + alfa) + Identity: a sin x + b cos x = sqrt(a**2 + b**2) sin (x + alpha) / | atan(b/a) for a > 0 - alfa < + alpha < | pi + atan(b/a) for a < 0 \ - => sin (x + alfa) = c / sqrt(a**2 + b**2) + => sin (x + alpha) = c / sqrt(a**2 + b**2) must use both branches for sin(x) = n: x = u (+ 2k*pi) | x = pi - u (+ 2k*pi) @@ -1809,23 +1824,20 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): ctht = np.cos(tht0) stht = np.sin(tht0) - nchi = np.c_[0., cchi, schi].T + nchi = np.c_[0.0, cchi, schi].T - gHat0_l = -np.vstack([ceta * ctht, - seta * ctht, - stht]) + gHat0_l = -np.vstack([ceta * ctht, seta * ctht, stht]) - a = cchi*ceta*ctht - b = -cchi*stht - c = stht + schi*seta*ctht + a = cchi * ceta * ctht + b = -cchi * stht + c = stht + schi * seta * ctht # form solution - abMag = np.sqrt(a*a + b*b) - assert np.all(abMag > 0), \ - "Beam vector specification is infeasible!" + abMag = np.sqrt(a * a + b * b) + assert np.all(abMag > 0), "Beam vector specification is infeasible!" phaseAng = np.arctan2(b, a) rhs = c / abMag - rhs[abs(rhs) > 1.] = np.nan + rhs[abs(rhs) > 1.0] = np.nan rhsAng = np.arcsin(rhs) # write ome angle output arrays (NaNs persist here) @@ -1849,9 +1861,7 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): tmp_eta = np.empty(numGood) tmp_gvec = gHat0_l[:, goodOnes] for i in range(numGood): - rchi = rotMatOfExpMap( - np.tile(ome_min[goodOnes][i], (3, 1)) * nchi - ) + rchi = rotMatOfExpMap(np.tile(ome_min[goodOnes][i], (3, 1)) * nchi) gHat_l = np.dot(rchi, tmp_gvec[:, i].reshape(3, 1)) tmp_eta[i] = np.arctan2(gHat_l[1], gHat_l[0]) eta_min[goodOnes] = tmp_eta @@ -1860,129 +1870,78 @@ def getFriedelPair(tth0, eta0, *ome0, **kwargs): # - ome1 is in RADIANS here # - convert and put into [-180, 180] ome1 = mapAngle( - mapAngle( - r2d*ome_min, [-180, 180], units='degrees' - ) + c1*ome0, [-180, 180], units='degrees' + mapAngle(r2d * ome_min, [-180, 180], units='degrees') + c1 * ome0, + [-180, 180], + units='degrees', ) # put eta1 in [-180, 180] - eta1 = mapAngle(r2d*eta_min, [-180, 180], units='degrees') + eta1 = mapAngle(r2d * eta_min, [-180, 180], units='degrees') if not outputDegrees: - ome1 = d2r * ome1 - eta1 = d2r * eta1 + ome1 *= d2r + eta1 *= d2r return ome1, eta1 -def getDparms(lp, lpTag, radians=True): +def getDparms( + lp: np.ndarray, lpTag: str, radians: Optional[bool] = True +) -> np.ndarray: """ Utility routine for getting dparms, that is the lattice parameters without symmetry -- 'triclinic' - """ - latVecOps = latticeVectors(lp, tag=lpTag, radians=radians) - return latVecOps['dparms'] - -def LoadFormFactorData(): - """ - Script to read in a csv file containing information relating the - magnitude of Q (sin(th)/lambda) to atomic form factor - - - Notes: - Atomic form factor data gathered from the International Tables of - Crystallography: - - P. J. Brown, A. G. Fox, E. N. Maslen, M. A. O'Keefec and B. T. M. Willis, - "Chapter 6.1. Intensity of diffracted intensities", International Tables - for Crystallography (2006). Vol. C, ch. 6.1, pp. 554-595 - """ - - dir1 = os.path.split(valunits.__file__) - dataloc = os.path.join(dir1[0], 'data', 'FormFactorVsQ.csv') - - data = np.zeros((62, 99), float) - - # FIXME: marked broken by DP - jj = 0 - with open(dataloc, 'rU') as csvfile: - datareader = csv.reader(csvfile, dialect=csv.excel) - for row in datareader: - ii = 0 - for val in row: - data[jj, ii] = float(val) - ii += 1 - jj += 1 - - sinThOverLamdaList = data[:, 0] - ffDataList = data[:, 1:] - - return sinThOverLamdaList, ffDataList - - -def RetrieveAtomicFormFactor(elecNum, magG, sinThOverLamdaList, ffDataList): - """ Interpolates between tabulated data to find the atomic form factor - for an atom with elecNum electrons for a given magnitude of Q - - USAGE: - - ff = RetrieveAtomicFormFactor(elecNum,magG,sinThOverLamdaList,ffDataList) - - INPUTS: - - 1) elecNum, (1 x 1 float) number of electrons for atom of interest - - 2) magG (1 x 1 float) magnitude of G - - 3) sinThOverLamdaList (n x 1 float ndarray) form factor data is tabulated - in terms of sin(theta)/lambda (A^-1). - - 3) ffDataList (n x m float ndarray) form factor data is tabulated in terms - of sin(theta)/lambda (A^-1). Each column corresponds to a different - number of electrons - - OUTPUTS: - - 1) ff (n x 1 float) atomic form factor for atom and hkl of interest - - NOTES: - Data should be calculated in terms of G at some point + Parameters + ---------- + lp : np.ndarray + Parsed lattice parameters + lpTag : str + Tag for the symmetry group of the lattice (from Laue group) + radians : bool, optional + Whether or not to use radians for angles, default is True + Returns + ------- + np.ndarray + The lattice parameters without symmetry. """ - sinThOverLambda = 0.5*magG - # lambda=2*d*sin(th) - # lambda=2*sin(th)/G - # 1/2*G=sin(th)/lambda - - ff = np.interp( - sinThOverLambda, sinThOverLamdaList, ffDataList[:, (elecNum - 1)] - ) - - return ff + latVecOps = latticeVectors(lp, tag=lpTag, radians=radians) + return latVecOps['dparms'] -def lorentz_factor(tth): +def lorentz_factor(tth: np.ndarray) -> np.ndarray: """ 05/26/2022 SS adding lorentz factor computation to the detector so that it can be compenstated for in the intensity correction - parameters: tth two theta of every pixel in radians + Parameters + ---------- + tth: np.ndarray + 2-theta of every pixel in radians + + Returns + ------- + np.ndarray + Lorentz factor for each pixel """ - theta = 0.5*tth + theta = 0.5 * tth cth = np.cos(theta) - sth2 = np.sin(theta)**2 + sth2 = np.sin(theta) ** 2 - L = 1./(4.0*cth*sth2) + return 1.0 / (4.0 * cth * sth2) - return L - -def polarization_factor(tth, unpolarized=True, eta=None, - f_hor=None, f_vert=None): +def polarization_factor( + tth: np.ndarray, + unpolarized: Optional[bool] = True, + eta: Optional[np.ndarray] = None, + f_hor: Optional[float] = None, + f_vert: Optional[float] = None, +) -> np.ndarray: """ 06/14/2021 SS adding lorentz polarization factor computation to the detector so that it can be compenstated for in the @@ -1998,13 +1957,16 @@ def polarization_factor(tth, unpolarized=True, eta=None, f_vert fraction of vertical polarization (~0 for XFELs) notice f_hor + f_vert = 1 + + FIXME, called without parameters like eta, f_hor, f_vert, but they default + to none in the current implementation, which will throw an error. """ - ctth2 = np.cos(tth)**2 + ctth2 = np.cos(tth) ** 2 if unpolarized: return (1 + ctth2) / 2 - seta2 = np.sin(eta)**2 - ceta2 = np.cos(eta)**2 - return f_hor*(seta2 + ceta2*ctth2) + f_vert*(ceta2 + seta2*ctth2) + seta2 = np.sin(eta) ** 2 + ceta2 = np.cos(eta) ** 2 + return f_hor * (seta2 + ceta2 * ctth2) + f_vert * (ceta2 + seta2 * ctth2) diff --git a/tests/planedata/test_exclusion.py b/tests/planedata/test_exclusion.py index bc8954b5c..f85c0fc47 100644 --- a/tests/planedata/test_exclusion.py +++ b/tests/planedata/test_exclusion.py @@ -24,10 +24,7 @@ def test_exclusion(): assert pd.nHKLs == 2 assert pd.getNhklRef() == 3 - try: - # Exclude a range - pd.exclusions = [[0, 2]] - assert pd.nHKLs == 1 - assert pd.getNhklRef() == 3 - except NotImplementedError: - pass + # Exclude a range + pd.exclusions = [[0, 2]] + assert pd.nHKLs == 1 + assert pd.getNhklRef() == 3 diff --git a/tests/planedata/test_misc.py b/tests/planedata/test_misc.py index 9a1686439..78a977039 100644 --- a/tests/planedata/test_misc.py +++ b/tests/planedata/test_misc.py @@ -1,3 +1,4 @@ +import os import numpy as np from hexrd.material.crystallography import PlaneData @@ -16,6 +17,8 @@ def test_misc(): 0.5, 0.01, ) + # Set ACK_DEPRECATED environment variable to true + os.environ['ACK_DEPRECATED'] = 'true' assert pd.laueGroup == 'C2h' assert pd.getLatticeType() == 'monoclinic' diff --git a/tests/planedata/test_with_data.py b/tests/planedata/test_with_data.py index 7b1e3ae61..b4e771792 100644 --- a/tests/planedata/test_with_data.py +++ b/tests/planedata/test_with_data.py @@ -1,3 +1,4 @@ +import os import numpy as np import pytest @@ -56,6 +57,7 @@ def flatten(xs): def test_plane_data_with_data(test_data_dir, materials): data = np.load(test_data_dir / 'plane_data_test.npy', allow_pickle=True) + os.environ['ACK_DEPRECATED'] = 'true' for obj in data: pd = materials[obj['mat_name']] assertEqualNumpyArrays(pd.hkls, obj['hkls']) @@ -84,6 +86,7 @@ def test_plane_data_with_data(test_data_dir, materials): # def test_with_data_maker(test_data_dir, materials): # data = [] +# os.environ['ACK_DEPRECATED'] = 'true' # for mat_name, pd in materials.items(): # obj = {"mat_name": mat_name} # obj['hkls'] = pd.hkls From 503d5f7d3d5f3cd0329d2222973ad4a574fe8ad2 Mon Sep 17 00:00:00 2001 From: Kevin Lewis <55998034+kevindlewis23@users.noreply.github.com> Date: Thu, 15 Aug 2024 10:35:48 -0400 Subject: [PATCH 02/16] Replace validateAngleRanges with xfcapi's version, should actually work on cw inputs now (#703) --- hexrd/xrdutil/utils.py | 104 +---------------------------------------- 1 file changed, 1 insertion(+), 103 deletions(-) diff --git a/hexrd/xrdutil/utils.py b/hexrd/xrdutil/utils.py index 1e9961dbb..2245dbdbc 100644 --- a/hexrd/xrdutil/utils.py +++ b/hexrd/xrdutil/utils.py @@ -262,109 +262,7 @@ def make_polar_net( return np.vstack(pts) -def validateAngleRanges( - angList: Union[np.ndarray, list[float]], - startAngs: Union[np.ndarray, list[float]], - stopAngs: Union[np.ndarray, list[float]], - ccw: bool = True, -) -> np.ndarray: - """ - Indetify angles that fall within specified ranges. - - 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 - """ - angList = np.atleast_1d(angList).flatten() - startAngs = np.atleast_1d(startAngs).flatten() - stopAngs = np.atleast_1d(stopAngs).flatten() - - if len(startAngs) != len(stopAngs): - raise ValueError("start and stop angles must have same length") - - # to avoid warnings in >=, <= later down, mark nans; - # need these to trick output to False in the case of nan input - nan_mask = np.isnan(angList) - - reflInRange = np.zeros(angList.shape, dtype=bool) - - # bin length for chunking - binLen = np.pi / 2.0 - - # in plane vectors defining wedges - x0 = np.vstack([np.cos(startAngs), np.sin(startAngs)]) - x1 = np.vstack([np.cos(stopAngs), np.sin(stopAngs)]) - - # dot products - dp = np.sum(x0 * x1, axis=0) - if np.any(dp >= 1.0 - sqrt_epsf) and len(startAngs) > 1: - # ambiguous case - raise RuntimeError( - "Improper usage; " - + "at least one of your ranges is already 360 degrees!" - ) - elif dp[0] >= 1.0 - sqrt_epsf and len(startAngs) == 1: - # trivial case! - reflInRange = np.ones(angList.shape, dtype=bool) - reflInRange[nan_mask] = False - else: - # solve for arc lengths - # ...note: no zeros should have made it here - a = x0[0, :] * x1[1, :] - x0[1, :] * x1[0, :] - b = x0[0, :] * x1[0, :] + x0[1, :] * x1[1, :] - phi = np.arctan2(b, a) - - arclen = 0.5 * np.pi - phi # these are clockwise - cw_phis = arclen < 0 - arclen[cw_phis] += 2 * np.pi # all positive (CW) now - if not ccw: - arclen = 2 * np.pi - arclen - - if sum(arclen) > 2 * np.pi: - raise RuntimeWarning( - "Specified angle ranges sum to > 360 degrees, " - + "which is suspect..." - ) - - # check that there are no more than dp = np.zeros(len(startAngs)) - for i in range(len(startAngs)): - # number or subranges using 'binLen' - numSubranges = int(np.ceil(arclen[i] / binLen)) - - # check remainder - binrem = np.remainder(arclen[i], binLen) - finalBinLen = binLen if binrem == 0 else binrem - - # if clockwise, negate bin length - if not ccw: - binLen = -binLen - finalBinLen = -finalBinLen - - # Create sub ranges on the fly to avoid ambiguity in dot product - # for wedges >= 180 degrees - subRanges = np.array( - [startAngs[i] + binLen * j for j in range(numSubranges)] - + [startAngs[i] + binLen * (numSubranges - 1) + finalBinLen] - ) - - for k in range(numSubranges): - zStart = _zproject(angList, subRanges[k]) - zStop = _zproject(angList, subRanges[k + 1]) - if ccw: - zStart[nan_mask] = 999.0 - zStop[nan_mask] = -999.0 - reflInRange = reflInRange | np.logical_and( - zStart <= 0, zStop >= 0 - ) - else: - zStart[nan_mask] = -999.0 - zStop[nan_mask] = 999.0 - reflInRange = reflInRange | np.logical_and( - zStart >= 0, zStop <= 0 - ) - return reflInRange +validateAngleRanges = xfcapi.validate_angle_ranges @deprecated(removal_date='2025-01-01') From 433a0d59f8da58c66192acd1b0bd0dea2bf9f5c1 Mon Sep 17 00:00:00 2001 From: Kevin Lewis <55998034+kevindlewis23@users.noreply.github.com> Date: Thu, 15 Aug 2024 10:36:39 -0400 Subject: [PATCH 03/16] Unitcell testcases (#702) * Unitcell testcases --- hexrd/material/unitcell.py | 177 ++++++++++++++++---------------- tests/unitcell/test_vec_math.py | 103 +++++++++++++++++++ 2 files changed, 194 insertions(+), 86 deletions(-) create mode 100644 tests/unitcell/test_vec_math.py diff --git a/hexrd/material/unitcell.py b/hexrd/material/unitcell.py index 9f3b3f243..26f3a2fef 100644 --- a/hexrd/material/unitcell.py +++ b/hexrd/material/unitcell.py @@ -42,7 +42,7 @@ def _calcstar(v, sym, mat): if dist < 1E-3: isnew = False break - if(isnew): + if isnew: vp = np.atleast_2d(vp) vsym = np.vstack((vsym, vp)) @@ -158,7 +158,7 @@ def calcmatrices(self): [a*c*cb, b*c*ca, c**2]]) self._vol = np.sqrt(np.linalg.det(self.dmt)) - if(self.vol < 1e-5): + if self.vol < 1e-5: warnings.warn('unitcell volume is suspiciously small') ''' @@ -198,32 +198,34 @@ def calcmatrices(self): choices are 'd' (direct), 'r' (reciprocal) and 'c' (cartesian)''' def TransSpace(self, v_in, inspace, outspace): - if(inspace == 'd'): - if(outspace == 'r'): + if inspace == outspace: + return v_in + if inspace == 'd': + if outspace == 'r': v_out = np.dot(v_in, self.dmt) - elif(outspace == 'c'): + elif outspace == 'c': v_out = np.dot(self.dsm, v_in) else: raise ValueError( - 'inspace in ''d'' but outspace can''t be identified') + 'inspace in "d" but outspace can\'t be identified') - elif(inspace == 'r'): - if(outspace == 'd'): + elif inspace == 'r': + if outspace == 'd': v_out = np.dot(v_in, self.rmt) - elif(outspace == 'c'): + elif outspace == 'c': v_out = np.dot(self.rsm, v_in) else: raise ValueError( - 'inspace in ''r'' but outspace can''t be identified') + 'inspace in "r" but outspace can\'t be identified') - elif(inspace == 'c'): - if(outspace == 'r'): - v_out = np.dot(v_in, self.rsm) - elif(outspace == 'd'): + elif inspace == 'c': + if outspace == 'r': v_out = np.dot(v_in, self.dsm) + elif outspace == 'd': + v_out = np.dot(v_in, self.rsm) else: raise ValueError( - 'inspace in ''c'' but outspace can''t be identified') + 'inspace in "c" but outspace can\'t be identified') else: raise ValueError('incorrect inspace argument') @@ -234,11 +236,11 @@ def TransSpace(self, v_in, inspace, outspace): def CalcDot(self, u, v, space): - if(space == 'd'): + if space == 'd': dot = np.dot(u, np.dot(self.dmt, v)) - elif(space == 'r'): + elif space == 'r': dot = np.dot(u, np.dot(self.rmt, v)) - elif(space == 'c'): + elif space == 'c': dot = np.dot(u, v) else: raise ValueError('space is unidentified') @@ -247,13 +249,13 @@ def CalcDot(self, u, v, space): def CalcLength(self, u, space): - if(space == 'd'): + if space == 'd': mat = self.dmt # vlen = np.sqrt(np.dot(u, np.dot(self.dmt, u))) - elif(space == 'r'): + elif space == 'r': mat = self.rmt # vlen = np.sqrt(np.dot(u, np.dot(self.rmt, u))) - elif(space == 'c'): + elif space == 'c': mat = np.eye(3) # vlen = np.linalg.norm(u) else: @@ -297,7 +299,7 @@ def CalcAngle(self, u, v, space): def CalcCross(self, p, q, inspace, outspace, vol_divide=False): iv = 0 - if(vol_divide): + if vol_divide: vol = self.vol else: vol = 1.0 @@ -306,50 +308,50 @@ def CalcCross(self, p, q, inspace, outspace, vol_divide=False): p[2]*q[0]-p[0]*q[2], p[0]*q[1]-p[1]*q[0]]) - if(inspace == 'd'): + if inspace == 'd': ''' cross product vector is in reciprocal space and can be converted to direct or cartesian space ''' pxq *= vol - if(outspace == 'r'): + if outspace == 'r': pass - elif(outspace == 'd'): + elif outspace == 'd': pxq = self.TransSpace(pxq, 'r', 'd') - elif(outspace == 'c'): + elif outspace == 'c': pxq = self.TransSpace(pxq, 'r', 'c') else: raise ValueError( 'inspace is ''d'' but outspace is unidentified') - elif(inspace == 'r'): + elif inspace == 'r': ''' cross product vector is in direct space and can be converted to any other space ''' pxq /= vol - if(outspace == 'r'): + if outspace == 'r': pxq = self.TransSpace(pxq, 'd', 'r') - elif(outspace == 'd'): + elif outspace == 'd': pass - elif(outspace == 'c'): + elif outspace == 'c': pxq = self.TransSpace(pxq, 'd', 'c') else: raise ValueError( 'inspace is ''r'' but outspace is unidentified') - elif(inspace == 'c'): + elif inspace == 'c': ''' cross product is already in cartesian space so no volume factor is involved. can be converted to any other space too ''' - if(outspace == 'r'): + if outspace == 'r': pxq = self.TransSpace(pxq, 'c', 'r') - elif(outspace == 'd'): + elif outspace == 'd': pxq = self.TransSpace(pxq, 'c', 'd') - elif(outspace == 'c'): + elif outspace == 'c': pass else: raise ValueError( @@ -398,7 +400,7 @@ def GenerateCartesianPGSym(self): self.SYM_PG_c = np.array(self.SYM_PG_c) self.SYM_PG_c[np.abs(self.SYM_PG_c) < eps] = 0. - if(self._pointGroup == self._laueGroup): + if self._pointGroup == self._laueGroup: self.SYM_PG_c_laue = self.SYM_PG_c else: for sop in self.SYM_PG_d_laue: @@ -422,8 +424,7 @@ def GenerateCartesianPGSym(self): supergroup_laue = self._supergroup_laue sym_supergroup_laue = symmetry.GeneratePGSYM(supergroup_laue) - if((self.latticeType == 'monoclinic' or - self.latticeType == 'triclinic')): + if self.latticeType in ('monoclinic', 'triclinic'): ''' for monoclinic groups c2 and c2h, the supergroups are orthorhombic, so no need to convert from direct to @@ -462,7 +463,7 @@ def GenerateCartesianPGSym(self): not be rotated as they have the c* axis already aligned with the z-axis SS 12/10/2020 ''' - if(self.latticeType == 'monoclinic'): + if self.latticeType == 'monoclinic': om = np.array([[1., 0., 0.], [0., 0., 1.], [0., -1., 0.]]) @@ -480,7 +481,7 @@ def GenerateCartesianPGSym(self): triclinic group c1!! SS 12/10/2020 ''' - if(self._pointGroup == 'c1'): + if self._pointGroup == 'c1': om = np.array([[1., 0., 0.], [0., 0., 1.], [0., -1., 0.]]) for i, s in enumerate(self.SYM_PG_supergroup): @@ -536,7 +537,7 @@ def CalcOrbit(self, v, reduceToUC=True): break # if its new add this to the list - if(isnew): + if isnew: asym_pos = np.vstack((asym_pos, rr)) n += 1 @@ -549,21 +550,21 @@ def CalcStar(self, v, space, applyLaue=False): this function calculates the symmetrically equivalent hkls (or uvws) for the reciprocal (or direct) point group symmetry. ''' - if(space == 'd'): + if space == 'd': mat = self.dmt.astype(np.float64) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_d_laue.astype(np.float64) else: sym = self.SYM_PG_d.astype(np.float64) - elif(space == 'r'): + elif space == 'r': mat = self.rmt.astype(np.float64) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_r_laue.astype(np.float64) else: sym = self.SYM_PG_r.astype(np.float64) - elif(space == 'c'): + elif space == 'c': mat = np.eye(3) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_c_laue.astype(np.float64) else: sym = self.SYM_PG_c.astype(np.float64) @@ -750,7 +751,11 @@ def InitializeInterpTable(self): f_anomalous_data = [] self.pe_cs = {} - data = importlib.resources.open_binary(hexrd.resources, 'Anomalous.h5') + data = ( + importlib.resources.files(hexrd.resources) + .joinpath('Anomalous.h5') + .open('rb') + ) with h5py.File(data, 'r') as fid: for i in range(0, self.atom_ntype): @@ -917,7 +922,7 @@ def ChooseSymmetric(self, hkllist, InversionSymmetry=True): laue = InversionSymmetry for i, g in enumerate(hkllist): - if(mask[i]): + if mask[i]: geqv = self.CalcStar(g, 'r', applyLaue=laue) @@ -995,11 +1000,11 @@ def getHKLs(self, dmin): for g in hkl_allowed: # ignore [0 0 0] as it is the direct beam - if(np.sum(np.abs(g)) != 0): + if np.sum(np.abs(g)) != 0: dspace = 1./self.CalcLength(g, 'r') - if(dspace >= dmin): + if dspace >= dmin: hkl_dsp.append(g) ''' @@ -1037,7 +1042,7 @@ def Required_C(self, C): return np.array([C[x] for x in _StiffnessDict[self._laueGroup][0]]) def MakeStiffnessMatrix(self, inp_Cvals): - if(len(inp_Cvals) != len(_StiffnessDict[self._laueGroup][0])): + if len(inp_Cvals) != len(_StiffnessDict[self._laueGroup][0]): x = len(_StiffnessDict[self._laueGroup][0]) msg = (f"number of constants entered is not correct." f" need a total of {x} independent constants.") @@ -1079,16 +1084,16 @@ def inside_spheretriangle(self, conn, dir3, hemisphere, switch): first get vertices of the triangles in the ''' vertex = self.sphere_sector.vertices[switch] - # if(switch == 'pg'): + # if switch == 'pg': # vertex = self.sphere_sector.vertices - # elif(switch == 'laue'): + # elif switch == 'laue': # vertex = self.sphere_sector.vertices_laue - # elif(switch == 'super'): + # elif switch == 'super': # vertex = self.sphere_sector.vertices_supergroup - # elif(switch == 'superlaue'): + # elif switch == 'superlaue': # vertex = self.sphere_sector.vertices_supergroup_laue A = np.atleast_2d(vertex[:, conn[0]]).T @@ -1107,29 +1112,29 @@ def inside_spheretriangle(self, conn, dir3, hemisphere, switch): determinant can be very small positive or negative number ''' - if(np.abs(d1) < eps): + if np.abs(d1) < eps: d1 = 0. - if(np.abs(d2) < eps): + if np.abs(d2) < eps: d2 = 0. - if(np.abs(d3) < eps): + if np.abs(d3) < eps: d3 = 0. ss = np.unique(np.sign([d1, d2, d3])) - if(hemisphere == 'upper'): - if(np.all(ss >= 0.)): + if hemisphere == 'upper': + if np.all(ss >= 0.): mask.append(True) else: mask.append(False) - elif(hemisphere == 'both'): - if(len(ss) == 1): + elif hemisphere == 'both': + if len(ss) == 1: mask.append(True) - elif(len(ss) == 2): - if(0 in ss): + elif len(ss) == 2: + if 0 in ss: mask.append(True) else: mask.append(False) - elif(len(ss) == 3): + elif len(ss) == 3: mask.append(False) mask = np.array(mask) @@ -1161,7 +1166,7 @@ def reduce_dirvector(self, dir3, switch='pg'): ''' idx = np.arange(dir3.shape[0], dtype=np.int32) dir3 = np.ascontiguousarray(np.atleast_2d(dir3)) - if(dir3.ndim != 2): + if dir3.ndim != 2: raise RuntimeError("reduce_dirvector: invalid shape of dir3 array") ''' @@ -1173,10 +1178,10 @@ def reduce_dirvector(self, dir3, switch='pg'): ''' eps = constants.sqrt_epsf - if(np.all(np.abs(np.linalg.norm(dir3, axis=1) - 1.0) < eps)): + if np.all(np.abs(np.linalg.norm(dir3, axis=1) - 1.0) < eps): dir3n = dir3 else: - if(np.all(np.linalg.norm(dir3) > eps)): + if np.all(np.linalg.norm(dir3) > eps): dir3n = dir3/np.tile(np.linalg.norm(dir3, axis=1), [3, 1]).T else: raise RuntimeError( @@ -1200,30 +1205,30 @@ def reduce_dirvector(self, dir3, switch='pg'): ntriangle = self.sphere_sector.ntriangle[switch] connectivity = self.sphere_sector.connectivity[switch] - if(switch == 'pg'): + if switch == 'pg': sym = self.SYM_PG_c - elif(switch == 'super'): + elif switch == 'super': sym = self.SYM_PG_supergroup - elif(switch == 'laue'): + elif switch == 'laue': sym = self.SYM_PG_c_laue - elif(switch == 'superlaue'): + elif switch == 'superlaue': sym = self.SYM_PG_supergroup_laue for sop in sym: - if(dir3_copy.size != 0): + if dir3_copy.size != 0: dir3_sym = np.dot(sop, dir3_copy.T).T mask = np.zeros(dir3_sym.shape[0]).astype(bool) - if(ntriangle == 0): - if(hemisphere == 'both'): + if ntriangle == 0: + if hemisphere == 'both': mask = np.ones(dir3_sym.shape[0], dtype=bool) - elif(hemisphere == 'upper'): + elif hemisphere == 'upper': mask = dir3_sym[:, 2] >= 0. else: for ii in range(ntriangle): @@ -1232,8 +1237,8 @@ def reduce_dirvector(self, dir3, switch='pg'): hemisphere, switch) mask = np.logical_or(mask, tmpmask) - if(np.sum(mask) > 0): - if(dir3_reduced.size != 0): + if np.sum(mask) > 0: + if dir3_reduced.size != 0: dir3_reduced = np.vstack( (dir3_reduced, dir3_sym[mask, :])) idx_red = np.hstack((idx_red, idx[mask])) @@ -1269,7 +1274,7 @@ class which correctly color the orientations for this crystal class. the replaceing barycenter to pi-theta) ''' - if(laueswitch == True): + if laueswitch: ''' this is the case where we color orientations based on the laue group of the crystal. this is always going to be the case with x-ray which @@ -1280,7 +1285,7 @@ class which correctly color the orientations for this crystal class. the dir3, switch='superlaue') switch = 'superlaue' - elif(laueswitch == False): + else: ''' follow the logic in the function description ''' @@ -1317,7 +1322,7 @@ def color_orientations(self, ''' first make sure that the rotation matric is size nx3x3 ''' - if(rmats.ndim == 2): + if rmats.ndim == 2: rmats = np.atleast_3d(rmats).T else: assert rmats.ndim == 3, "rotations matrices need to \ @@ -1358,7 +1363,7 @@ def convert_lp_to_valunits(self, lp): """ lp_valunit = [] for i in range(6): - if(i < 3): + if i < 3: lp_valunit.append( valWUnit('lp', 'length', lp[i], 'nm')) @@ -1415,7 +1420,7 @@ def lparms(self, lp): self.calcmatrices() self.init_max_g_index() self.CalcMaxGIndex() - if(hasattr(self, 'numat')): + if hasattr(self, 'numat'): self.CalcDensity() @property @@ -1541,7 +1546,7 @@ def U(self): def U(self, Uarr): self._U = Uarr self.aniU = False - if(Uarr.ndim > 1): + if Uarr.ndim > 1: self.aniU = True self.calcBetaij() @@ -1569,9 +1574,9 @@ def sgnum(self): @sgnum.setter def sgnum(self, val): - if(not(isinstance(val, int))): + if not(isinstance(val, int)): raise ValueError('space group should be integer') - if(not((val >= 1) and (val <= 230))): + if not((val >= 1) and (val <= 230)): raise ValueError('space group number should be between 1 and 230.') self._sym_sgnum = val diff --git a/tests/unitcell/test_vec_math.py b/tests/unitcell/test_vec_math.py new file mode 100644 index 000000000..794ad4da5 --- /dev/null +++ b/tests/unitcell/test_vec_math.py @@ -0,0 +1,103 @@ +from pytest import fixture +from hexrd.material import Material, unitcell +import numpy as np + + +@fixture +def cell() -> unitcell.unitcell: + return Material().unitcell + + +def get_space(): + return np.random.choice(['d', 'r', 'c']) + + +def test_calc_dot(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + v2c = cell.TransSpace(v2, space, 'c') + assert np.allclose(np.dot(v1c, v2c), cell.CalcDot(v1, v2, space)) + + +def test_calc_length(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + assert np.allclose(np.linalg.norm(v1c), cell.CalcLength(v1, space)) + + +def test_calc_angle(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + v2c = cell.TransSpace(v2, space, 'c') + norms = np.linalg.norm(v1c) * np.linalg.norm(v2c) + assert np.allclose( + np.arccos(np.dot(v1c, v2c) / norms), cell.CalcAngle(v1, v2, space) + ) + + +def test_calc_cross(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + inspace = get_space() + outspace = get_space() + v1c = cell.TransSpace(v1, inspace, 'c') + v2c = cell.TransSpace(v2, inspace, 'c') + expected = cell.TransSpace(np.cross(v1c, v2c), 'c', outspace) + assert np.allclose( + expected, cell.CalcCross(v1, v2, inspace, outspace, True) + ) + + +def test_norm_vec(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = get_space() + norm_v1 = cell.NormVec(v1, space) + assert np.allclose(1, cell.CalcLength(norm_v1, space)) + # Make sure we don't change the direction + assert np.allclose( + v1 / np.linalg.norm(v1), norm_v1 / np.linalg.norm(norm_v1) + ) + + +def test_trans_space(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + inspace = get_space() + outspace = get_space() + v1_out = cell.TransSpace(v1, inspace, outspace) + assert np.allclose(v1, cell.TransSpace(v1_out, outspace, inspace)) + + +def test_calc_star(cell: unitcell.unitcell): + """ + Just ensuring that the outspace doesn't matter + """ + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = np.random.choice(['d', 'r']) + v1c = cell.TransSpace(v1, space, 'c') + assert np.allclose( + cell.CalcStar(v1, space, False), + cell.TransSpace(cell.CalcStar(v1c, 'c', False), 'c', space), + ) + assert np.allclose( + cell.CalcStar(v1, space, True), + cell.TransSpace(cell.CalcStar(v1c, 'c', True), 'c', space), + ) From 5a7f680da1ed0ea9c7483c472324ca4805de3f0c Mon Sep 17 00:00:00 2001 From: Kevin Lewis Date: Thu, 15 Aug 2024 14:22:14 -0400 Subject: [PATCH 04/16] Calibration testcases --- hexrd/matrixutil.py | 1 + tests/data/calibration_expected.npy | Bin 0 -> 940 bytes tests/test_calibration.py | 183 ++++++++++++++++++++++++++++ 3 files changed, 184 insertions(+) create mode 100644 tests/data/calibration_expected.npy create mode 100644 tests/test_calibration.py diff --git a/hexrd/matrixutil.py b/hexrd/matrixutil.py index 64ce134ed..d48c9f87d 100644 --- a/hexrd/matrixutil.py +++ b/hexrd/matrixutil.py @@ -182,6 +182,7 @@ def vecMVToSymm(A, scale=True): convert from Mandel-Voigt vector to symmetric matrix representation (JVB) """ + A = np.atleast_1d(A).flatten() if scale: fac = sqr2 else: diff --git a/tests/data/calibration_expected.npy b/tests/data/calibration_expected.npy new file mode 100644 index 0000000000000000000000000000000000000000..e4e0aa5bbcd3c992cc59257e26a74b287f498ae4 GIT binary patch literal 940 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1Jt-^H=`z0QA$Z=K`K`vTLcpW1B1UsA$w;> zdm%?qA*Y5na|9z$tfr95&(F{6KM;TkZ~Kx$?xfDxLY~?}UX2Jgppv58#FBWULcR=! z46s%F89;jkf(iwJ_INXUGq)8AWk4)QvnUh>Ss)TrC26G)w6YmWFr zH~R;ZC60Xl|9*c%q4$pX^z8kG5+#L_NeGvk6iUHdDh+a}Oi-aL+@*3EqCl%xZD4+r z^Y{IJ2ak76_FU%s54@Ng!?1eh{SR>$wx0*94J&PDO5)?$fQsMW{@VxAg!Q6ZKy#yV1p{> zFWgia^4R`(kHx(b2Mvc?59EZHCjt%9Eh*GXir@rBTS{hPZhl@$d}5(Kj&xO=oS2hX z#8qgJA?nQG!xOWZvx&!{(6Fu0D5%glf(K|udQoC#UVK4fQDSa!p$W`MrXVMo1r?fu zoW#>wXaRCkLiSwA(l=Kf4%~fcU#740!Cw0AoYra4J`O-*@PH5Y5J_|zQ%j*`NugCz LX<|vCb&?(cl!G!k literal 0 HcmV?d00001 diff --git a/tests/test_calibration.py b/tests/test_calibration.py new file mode 100644 index 000000000..600d53983 --- /dev/null +++ b/tests/test_calibration.py @@ -0,0 +1,183 @@ +import h5py +import numpy as np +import yaml + +import pytest + +from hexrd.material.material import load_materials_hdf5 +from hexrd.instrument.hedm_instrument import HEDMInstrument +from hexrd import rotations as rot + +from hexrd.fitting.calibration import ( + InstrumentCalibrator, + LaueCalibrator, + PowderCalibrator, +) + + +@pytest.fixture +def calibration_dir(example_repo_path): + return example_repo_path / 'tardis' / 'simulated' + + +def test_calibration(calibration_dir, test_data_dir): + # Load the material + with h5py.File(calibration_dir / 'materials.h5', 'r') as rf: + materials = load_materials_hdf5(rf) + + # Load the picks + with open( + calibration_dir / 'uncalibrated_tardis.yml', 'r', encoding='utf-8' + ) as rf: + conf = yaml.safe_load(rf) + + instrument = HEDMInstrument(conf) + + # Load the images + img_npz = np.load(calibration_dir / 'tardis_images.npz') + img_dict = {k: img_npz[k] for k in img_npz} + + # Load picks + picks = np.load(calibration_dir / 'picks.npy', allow_pickle=True) + + euler_convention = {'axes_order': 'zxz', 'extrinsic': False} + + # Create the calibrators + calibrators = [] + for pick_data in picks: + if pick_data['type'] == 'powder': + kwargs = { + 'instr': instrument, + 'material': materials[pick_data['material']], + 'img_dict': img_dict, + 'default_refinements': pick_data['default_refinements'], + 'tth_distortion': pick_data['tth_distortion'], + 'calibration_picks': pick_data['picks'], + } + calibrators.append(PowderCalibrator(**kwargs)) + + elif pick_data['type'] == 'laue': + # gpflags = [i[1] for i in pick_data['refinements']] + kwargs = { + 'instr': instrument, + 'material': materials[pick_data['material']], + 'grain_params': pick_data['options']['crystal_params'], + 'default_refinements': pick_data['default_refinements'], + 'min_energy': pick_data['options']['min_energy'], + 'max_energy': pick_data['options']['max_energy'], + 'calibration_picks': pick_data['picks'], + 'euler_convention': euler_convention, + } + calibrators.append(LaueCalibrator(**kwargs)) + + calibrator = InstrumentCalibrator( + *calibrators, + engineering_constraints='TARDIS', + euler_convention=euler_convention, + ) + + x0 = calibrator.params.valuesdict() + result = calibrator.run_calibration(None) + x1 = result.params.valuesdict() + + # Parse the data + tilt_angle_names = [ + [ + f'IMAGE_PLATE_{n}_euler_z', + f'IMAGE_PLATE_{n}_euler_xp', + f'IMAGE_PLATE_{n}_euler_zpp', + ] + for n in [2, 4] + ] + + rmats = { + 'old': [ + euler_to_rot([x0[k] for k in names]) for names in tilt_angle_names + ], + 'new': [ + euler_to_rot([x1[k] for k in names]) for names in tilt_angle_names + ], + } + + tvec_names = [ + [ + f'IMAGE_PLATE_{n}_tvec_x', + f'IMAGE_PLATE_{n}_tvec_y', + f'IMAGE_PLATE_{n}_tvec_z', + ] + for n in [2, 4] + ] + + tvecs = { + 'old': [np.array([x0[k] for k in vec_names]) for vec_names in tvec_names], + 'new': [np.array([x1[k] for k in vec_names]) for vec_names in tvec_names], + } + + grain_param_names = [f'LiF_grain_param_{n}' for n in range(12)] + grain_params = { + 'old': np.array([x0[name] for name in grain_param_names]), + 'new': np.array([x1[name] for name in grain_param_names]), + } + + diamond_a_vals = { + 'old': x0['diamond_a'], + 'new': x1['diamond_a'], + } + + expected = np.load( + test_data_dir / 'calibration_expected.npy', allow_pickle=True + ) + + assert_errors_are_better( + rmats, tvecs, grain_params, diamond_a_vals, expected.item() + ) + + +def euler_to_rot(euler): + return rot.RotMatEuler(np.array(euler), 'zxz', False, 'degrees').rmat + + +def assert_errors_are_better( + rmats, tvecs, grain_params, diamond_a_vals, expected +): + """ + Make sure error has decreased during fitting + """ + # What fraction of the old error we need to have (at worst) for the + # test to pass + ERROR_TOLERANCE = 0.5 + + rmat_error_2_old = np.linalg.norm( + rmats['old'][0] @ expected['rmat_2'].T - np.eye(3) + ) + rmat_error_2_new = np.linalg.norm( + rmats['new'][0] @ expected['rmat_2'].T - np.eye(3) + ) + rmat_error_4_old = np.linalg.norm( + rmats['old'][1] @ expected['rmat_4'].T - np.eye(3) + ) + rmat_error_4_new = np.linalg.norm( + rmats['new'][1] @ expected['rmat_4'].T - np.eye(3) + ) + assert rmat_error_2_new < rmat_error_2_old * ERROR_TOLERANCE + assert rmat_error_4_new < rmat_error_4_old * ERROR_TOLERANCE + + tvec_error_2_old = np.linalg.norm(tvecs['old'][0] - expected['tvec_2']) + tvec_error_2_new = np.linalg.norm(tvecs['new'][0] - expected['tvec_2']) + tvec_error_4_old = np.linalg.norm(tvecs['old'][1] - expected['tvec_4']) + tvec_error_4_new = np.linalg.norm(tvecs['new'][1] - expected['tvec_4']) + + assert tvec_error_2_new < tvec_error_2_old * ERROR_TOLERANCE + assert tvec_error_4_new < tvec_error_4_old * ERROR_TOLERANCE + + grain_param_error_old = np.linalg.norm( + grain_params['old'] - expected['grain_params'] + ) + grain_param_error_new = np.linalg.norm( + grain_params['new'] - expected['grain_params'] + ) + assert grain_param_error_new < grain_param_error_old * ERROR_TOLERANCE + + diamond_a_error_old = np.abs(diamond_a_vals['old'] - expected['diamond_a']) + diamond_a_error_new = np.abs(diamond_a_vals['new'] - expected['diamond_a']) + assert diamond_a_error_new < diamond_a_error_old * ERROR_TOLERANCE From 58fd4e09cdce12daf692dc2ce3da39e0946a8571 Mon Sep 17 00:00:00 2001 From: Kevin Lewis Date: Thu, 15 Aug 2024 15:15:17 -0400 Subject: [PATCH 05/16] More accurate calibration. --- tests/data/calibration_expected.npy | Bin 940 -> 940 bytes tests/test_calibration.py | 15 +++++++++++++-- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/tests/data/calibration_expected.npy b/tests/data/calibration_expected.npy index e4e0aa5bbcd3c992cc59257e26a74b287f498ae4..5dadf135d816539702f6f201939623ab20977e20 100644 GIT binary patch delta 357 zcmV-r0h<1-2doFMZvjWijG*85`tLvZsHrGU32Hy+BrKU&-l0E~Uf84^-l4w$5CFof z;6uMIkMlMH`|m%9k9=jMF={`Qm!$>u|L?y5C`t|?J3y1R0W?Pr#8;wO{qMhEprMRK z8D&30GZOOB=$XG|B_hl5@|i#Jje1nlrjS4XstG7b{O`Yv#}vKuw_(3G;&H#@{_npp z>J>tbwvdxC0x}kkX>+Vm^u528O@MO|VRB(@b8#&@~dY}K$Bqu6q9fR1(R?CFcma( zaW!IC2><{9WNBe-Z*F8?VR1H-+5#UOH)s_@2zeJ7^#6_wKyf&RaXC_PI(1=maXVrz DGdG#` delta 348 zcmV-i0i*t`2doFMZvjU-C?=Z?`|m%%Gn}-t&O$#x7K;7A$)P`_u=B%t$)Uf1ZB09L zm21D)ZSOS#`|m%%A%$!gp+-O8lNiMF|L?zmaZ9*gXKRzT0W?RdumkLD{qMg(NbZ6^ z3OB#N=$%{utC_#e7B>1|sF^>fq(xyQ@sU5Y+o*{v{O`XYVnhT4=t)0JQ-kik{_nqe z_+zkR?U9o)0x}jGAdTM&p~An0D$)RV6eqygu_$B3(UW-sGcI|bp|W{Y<3GoYINo$X zC_via9ufnQIKXi&ba5_XSP1|C0Ay)lZEtR5Utw`ClW_wdT`*`BLkLU^T%!qs3_x)) zhH)}daWhy90001Ia$#w1UvOb^VQq79G-v>50ZRadaW!ZIhjBJiaW^PS0ZR;paX4rd ui9=v(of&oP)j+`A;XigSBJe*O-JFG%TTGJ@14;}zba6UjbzyXIJ7O+Hqn7Fb diff --git a/tests/test_calibration.py b/tests/test_calibration.py index 600d53983..1182f6ae2 100644 --- a/tests/test_calibration.py +++ b/tests/test_calibration.py @@ -77,7 +77,7 @@ def test_calibration(calibration_dir, test_data_dir): ) x0 = calibrator.params.valuesdict() - result = calibrator.run_calibration(None) + result = calibrator.run_calibration({'max_nfev': 1300}) x1 = result.params.valuesdict() # Parse the data @@ -124,6 +124,17 @@ def test_calibration(calibration_dir, test_data_dir): 'new': x1['diamond_a'], } + # expected_obj = { + # 'rmat_2': rmats['new'][0], + # 'rmat_4': rmats['new'][1], + # 'tvec_2': tvecs['new'][0], + # 'tvec_4': tvecs['new'][1], + # 'grain_params': grain_params['new'], + # 'diamond_a': diamond_a_vals['new'], + # } + + # np.save(test_data_dir / 'calibration_expected.npy', expected_obj) + expected = np.load( test_data_dir / 'calibration_expected.npy', allow_pickle=True ) @@ -145,7 +156,7 @@ def assert_errors_are_better( """ # What fraction of the old error we need to have (at worst) for the # test to pass - ERROR_TOLERANCE = 0.5 + ERROR_TOLERANCE = 0.8 rmat_error_2_old = np.linalg.norm( rmats['old'][0] @ expected['rmat_2'].T - np.eye(3) From 5145a0a32ab99ab9981d5c4d9ebbae3e643ebf80 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Thu, 15 Aug 2024 13:08:05 -0700 Subject: [PATCH 06/16] fixing some translation issues in origin choice 2 --- hexrd/material/symmetry.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/hexrd/material/symmetry.py b/hexrd/material/symmetry.py index ed8830fb3..5ea30a1bb 100644 --- a/hexrd/material/symmetry.py +++ b/hexrd/material/symmetry.py @@ -80,6 +80,7 @@ def GeneratorString(sgnum): def MakeGenerators(genstr, setting): + t = 'aOOO' mat = SYM_fillgen(t) genmat = mat @@ -113,22 +114,29 @@ def MakeGenerators(genstr, setting): if(genstr[istop] != '0'): if(setting != 0): t = genstr[istop+1:istop+4] - trans = np.array([constants.SYM_GENERATORS[t[0]],\ - constants.SYM_GENERATORS[t[1]],\ - constants.SYM_GENERATORS[t[2]] - ]) - for i in range(genmat.shape[0]): - genmat[i,0:3,3] -= trans + t = 'a' + t # get the translation without any rotation + sym = np.squeeze(SYM_fillgen(t, sgn=-1)) + sym2 = np.squeeze(SYM_fillgen(t)) + for i in range(1, genmat.shape[0]): + generator = np.dot(sym2, np.dot( + np.squeeze(genmat[i,:,:]), + sym)) + frac = np.modf(generator[0:3,3])[0] + frac[frac < 0.] += 1. + frac[np.abs(frac) < 1E-5] = 0.0 + frac[np.abs(frac-1.0) < 1E-5] = 0.0 + generator[0:3,3] = frac + genmat[i,:,:] = generator return genmat, centrosymmetric -def SYM_fillgen(t): +def SYM_fillgen(t, sgn=1): mat = np.zeros([4,4]) mat[3,3] = 1. mat[0:3,0:3] = constants.SYM_GENERATORS[t[0]] - mat[0:3,3] = np.array([constants.SYM_GENERATORS[t[1]],\ - constants.SYM_GENERATORS[t[2]],\ + mat[0:3,3] = sgn*np.array([constants.SYM_GENERATORS[t[1]], + constants.SYM_GENERATORS[t[2]], constants.SYM_GENERATORS[t[3]] ]) From 6a975de7fb0e43f8bddf88496bb02cd88acd42ca Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Thu, 15 Aug 2024 13:28:03 -0700 Subject: [PATCH 07/16] add space groups with 2 origin choices and the respective origin site symmetries. --- hexrd/constants.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/hexrd/constants.py b/hexrd/constants.py index aa813e52f..42cf1f0b2 100644 --- a/hexrd/constants.py +++ b/hexrd/constants.py @@ -1100,6 +1100,38 @@ def is_writable_file(path): SYM_GENERATORS['Y'] = -1./4. SYM_GENERATORS['Z'] = -1./8. +''' this dictionary contains the space group numbers +which have two origin choices the two values are the +site symmetries of the origin. There are 24 such space +groups +''' +two_origin_choice = { +48 : ['222', '-1'], +50 : ['222/n', '-1'], +59 : ['mm2/n', '-1'], +68 : ['222', '-1'], +70 : ['222', '-1'], +85 : ['-4', '-1'], +86 : ['-4', '-1'], +88 : ['-4', '-1'], +125: ['422', '2/m'], +126: ['422/n', '-1'], +129: ['-4m2', '2/m'], +130: ['-4/ncn', '-1'], +133: ['-4121/c', '-1'], +134: ['-42m', '2/m'], +137: ['-4m2/n', '-1'], +138: ['-4cg', '2/m'], +141: ['-4m2', '2/m'], +142: ['-4c21', '-1'], +201: ['23', '-3'], +203: ['23', '-3'], +222: ['432', '-3'], +224: ['-43m', '-3m'], +227: ['-43m', '-3m'], +228: ['23', '-3'], +} + ''' @AUTHOR Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov From 8655bd461199b1f50b8e60948b2c1413b4999f6f Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Thu, 15 Aug 2024 14:24:19 -0700 Subject: [PATCH 08/16] correctly handle passing of sgsetting for file reads --- hexrd/material/material.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/hexrd/material/material.py b/hexrd/material/material.py index 30b89e9b3..5ca1ca67a 100644 --- a/hexrd/material/material.py +++ b/hexrd/material/material.py @@ -37,7 +37,10 @@ from hexrd.material.crystallography import PlaneData as PData from hexrd.material import symmetry, unitcell from hexrd.valunits import valWUnit -from hexrd.constants import ptable, ptableinverse, chargestate +from hexrd.constants import (ptable, + ptableinverse, + chargestate, + two_origin_choice) from os import path from pathlib import Path @@ -76,6 +79,11 @@ def _kev(x): def _key(x): return x.name +def get_default_sgsetting(sgnum): + if sgnum in two_origin_choice: + return 1 + else: + return 0 # # ---------------------------------------------------CLASS: Material @@ -148,7 +156,7 @@ def __init__( material_file=None, dmin=None, kev=None, - sgsetting=DFLT_SGSETTING, + sgsetting=None, ): """Constructor for Material @@ -165,8 +173,6 @@ def __init__( self.description = '' - self.sgsetting = sgsetting - if material_file: # >> @ date 08/20/2020 SS removing dependence on hklmax if isinstance(material_file, (Path, str)): @@ -180,6 +186,9 @@ def __init__( self._readCif(material_file) elif isinstance(material_file, h5py.Group) or form in h5_suffixes: self._readHDFxtal(fhdf=material_file, xtal=name) + if sgsetting is not None: + if sgsetting in [0, 1]: + self._sgsetting = sgsetting else: # default name self._name = Material.DFLT_XTAL @@ -189,8 +198,8 @@ def __init__( # self.description = '' # - self.sgnum = Material.DFLT_SGNUM self._sgsetting = Material.DFLT_SGSETTING + self.sgnum = Material.DFLT_SGNUM # self._atominfo = Material.DFLT_ATOMINFO # @@ -645,6 +654,7 @@ def _readCif(self, fcif=DFLT_NAME + '.cif'): lparms[i] = _degrees(lparms[i]) self._lparms = lparms + self._sgsetting = get_default_sgsetting(sgnum) self.sgnum = sgnum # fractional atomic site, occ and vibration amplitude @@ -795,7 +805,6 @@ def _readCif(self, fcif=DFLT_NAME + '.cif'): self._atomtype = np.asarray(atomtype).astype(np.int32) self._charge = charge - self._sgsetting = 0 self._dmin = Material.DFLT_DMIN self._beamEnergy = Material.DFLT_KEV @@ -864,6 +873,9 @@ def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME): self._lparms = lparms # fill space group and lattice parameters + self._sgsetting = np.array( + gid.get('SpaceGroupSetting'), dtype=np.int32 + ).item() self.sgnum = sgnum # the U factors are related to B by the relation B = 8pi^2 U @@ -880,10 +892,6 @@ def _readHDFxtal(self, fhdf=DFLT_NAME, xtal=DFLT_NAME): self._charge = ['0'] * self._atomtype.shape[0] self._atom_ntype = self._atomtype.shape[0] - self._sgsetting = np.array( - gid.get('SpaceGroupSetting'), dtype=np.int32 - ).item() - if 'stiffness' in gid: # we're assuming the stiffness is in units of GPa self.stiffness = np.array(gid.get('stiffness')) @@ -1428,8 +1436,7 @@ def loadMaterialList(cfgFile): def load_materials_hdf5( - f, dmin=None, kev=None, sgsetting=Material.DFLT_SGSETTING -): + f, dmin=None, kev=None): """Load materials from an HDF5 file The file uses the HDF5 file format. @@ -1452,7 +1459,6 @@ def load_materials_hdf5( 'material_file': f, 'dmin': dmin, 'kev': kev, - 'sgsetting': sgsetting, } return {name: Material(name, **kwargs) for name in names} From 4118a8fecc9ff24d9fd626dbe9e61dd84b9828c9 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Thu, 15 Aug 2024 18:34:38 -0500 Subject: [PATCH 09/16] Fix up calibration test This fixes the expected output to match what the correct answers are for the simulated TARDIS. This also gets the lattice parameter to refine (it was already set to the ideal value before), and it removes the tilt refinements (we were refining too many things at once - this setup with just translations and materials seems to be working well). Signed-off-by: Patrick Avery --- tests/data/calibration_expected.npy | Bin 940 -> 938 bytes tests/test_calibration.py | 71 +++++++++------------------- 2 files changed, 23 insertions(+), 48 deletions(-) diff --git a/tests/data/calibration_expected.npy b/tests/data/calibration_expected.npy index 5dadf135d816539702f6f201939623ab20977e20..587793d8f5d7dec5448c05a7ea82330db4155527 100644 GIT binary patch delta 398 zcmZ3(zKVTAKBET%7<{m2ixDwaTeffyJA}4@@gedXY+&l3@_V5C{V@5-+Zl~vhV4gH z1X04kkYEW_4A%^^48t%ZCPQ(UHu=`qfjjQ09}rpR{;*2Ka&i@u5xSfMOztp~Nj+Nx z6HsfWMWIP&M|+`ZP@$QIH*d r#45BXDYQ(A-~=*LG81$2^HSmy3$5H8csiNxF6|9Wr#W5r4hKh7-s7b-e4BknDXvl-Lx_x2tAHK|Jsv+So#U&{LR|NH$68eV*g)(-n8Z)Y_2;5!n&DDv<7{Rs;e z_P9!>+B+EweY*5w=KfSwrPCij&9wj6TNQF?*#!Ilt2i}0|GeMdb6jNa=k1C6%^wx+ zfBf(Le*IUXPQBYE?4NAJWGLR3S-d9j%ijG{ecPnf1k?{CXe#ayyKgbMipfa4KkeHG zhS(?e?o;mAw>fAy$X_Vi=P=b`@?j^rPm;Ca-kH?|VqOH&}sL-l3v82#CNe=+Q3z`!E diff --git a/tests/test_calibration.py b/tests/test_calibration.py index 1182f6ae2..1939c96b4 100644 --- a/tests/test_calibration.py +++ b/tests/test_calibration.py @@ -1,3 +1,5 @@ +import copy + import h5py import numpy as np import yaml @@ -6,7 +8,6 @@ from hexrd.material.material import load_materials_hdf5 from hexrd.instrument.hedm_instrument import HEDMInstrument -from hexrd import rotations as rot from hexrd.fitting.calibration import ( InstrumentCalibrator, @@ -76,11 +77,6 @@ def test_calibration(calibration_dir, test_data_dir): euler_convention=euler_convention, ) - x0 = calibrator.params.valuesdict() - result = calibrator.run_calibration({'max_nfev': 1300}) - x1 = result.params.valuesdict() - - # Parse the data tilt_angle_names = [ [ f'IMAGE_PLATE_{n}_euler_z', @@ -89,16 +85,22 @@ def test_calibration(calibration_dir, test_data_dir): ] for n in [2, 4] ] + all_tilt_angle_names = tilt_angle_names[0] + tilt_angle_names[1] - rmats = { - 'old': [ - euler_to_rot([x0[k] for k in names]) for names in tilt_angle_names - ], - 'new': [ - euler_to_rot([x1[k] for k in names]) for names in tilt_angle_names - ], - } + # The tilts are already ideal. Do not refine. + orig_params = copy.deepcopy(calibrator.params) + for name in all_tilt_angle_names: + calibrator.params[name].vary = False + # The lattice parameter is actually perfect. Adjust it a little + # So that it can be corrected. + calibrator.params['diamond_a'].value = 3.58 + + x0 = calibrator.params.valuesdict() + result = calibrator.run_calibration({}) + x1 = result.params.valuesdict() + + # Parse the data tvec_names = [ [ f'IMAGE_PLATE_{n}_tvec_x', @@ -124,54 +126,24 @@ def test_calibration(calibration_dir, test_data_dir): 'new': x1['diamond_a'], } - # expected_obj = { - # 'rmat_2': rmats['new'][0], - # 'rmat_4': rmats['new'][1], - # 'tvec_2': tvecs['new'][0], - # 'tvec_4': tvecs['new'][1], - # 'grain_params': grain_params['new'], - # 'diamond_a': diamond_a_vals['new'], - # } - - # np.save(test_data_dir / 'calibration_expected.npy', expected_obj) - expected = np.load( test_data_dir / 'calibration_expected.npy', allow_pickle=True ) assert_errors_are_better( - rmats, tvecs, grain_params, diamond_a_vals, expected.item() + tvecs, grain_params, diamond_a_vals, expected.item() ) -def euler_to_rot(euler): - return rot.RotMatEuler(np.array(euler), 'zxz', False, 'degrees').rmat - - def assert_errors_are_better( - rmats, tvecs, grain_params, diamond_a_vals, expected + tvecs, grain_params, diamond_a_vals, expected ): """ Make sure error has decreased during fitting """ # What fraction of the old error we need to have (at worst) for the - # test to pass - ERROR_TOLERANCE = 0.8 - - rmat_error_2_old = np.linalg.norm( - rmats['old'][0] @ expected['rmat_2'].T - np.eye(3) - ) - rmat_error_2_new = np.linalg.norm( - rmats['new'][0] @ expected['rmat_2'].T - np.eye(3) - ) - rmat_error_4_old = np.linalg.norm( - rmats['old'][1] @ expected['rmat_4'].T - np.eye(3) - ) - rmat_error_4_new = np.linalg.norm( - rmats['new'][1] @ expected['rmat_4'].T - np.eye(3) - ) - assert rmat_error_2_new < rmat_error_2_old * ERROR_TOLERANCE - assert rmat_error_4_new < rmat_error_4_old * ERROR_TOLERANCE + # test to pass. For now, just make sure the error decreased. + ERROR_TOLERANCE = 1 tvec_error_2_old = np.linalg.norm(tvecs['old'][0] - expected['tvec_2']) tvec_error_2_new = np.linalg.norm(tvecs['new'][0] - expected['tvec_2']) @@ -191,4 +163,7 @@ def assert_errors_are_better( diamond_a_error_old = np.abs(diamond_a_vals['old'] - expected['diamond_a']) diamond_a_error_new = np.abs(diamond_a_vals['new'] - expected['diamond_a']) + + # The old diamond setting was actually perfect, but we let it refine + # The new diamond error should be less than 2% assert diamond_a_error_new < diamond_a_error_old * ERROR_TOLERANCE From 5602363152ff9e7ae91c313f26c37f160578a031 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Tue, 20 Aug 2024 11:37:36 -0700 Subject: [PATCH 10/16] PEP8 fixes. Moved two_origin_choice to materials.symbols --- hexrd/constants.py | 90 +++++++++++++------------------ hexrd/material/material.py | 13 ++--- hexrd/material/mksupport.py | 22 ++++---- hexrd/material/symbols.py | 42 +++++++++++---- hexrd/material/symmetry.py | 105 ++++++++++++++++++------------------ 5 files changed, 143 insertions(+), 129 deletions(-) diff --git a/hexrd/constants.py b/hexrd/constants.py index 42cf1f0b2..932814ae6 100644 --- a/hexrd/constants.py +++ b/hexrd/constants.py @@ -82,7 +82,8 @@ zeros_3x1 = np.zeros((3, 1)) zeros_6x1 = np.zeros((6, 1)) -# reference beam direction and eta=0 ref in LAB FRAME for standard geometry +'''reference beam direction and +eta=0 ref in LAB FRAME for standard geometry''' beam_vec = -lab_z eta_vec = lab_x @@ -162,7 +163,9 @@ 3628800.]).astype(np.complex128) """ ->> @AUTHOR: Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov +>> @AUTHOR: Saransh Singh, + Lawrence Livermore National Lab, + saransh1@llnl.gov >> @DATE: 11/28/2022 SS 1.0 original >> @DETAILS: constants for rodrigues FZ """ @@ -177,7 +180,9 @@ ) ''' ->> @AUTHOR: Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov +>> @AUTHOR: Saransh Singh, + Lawrence Livermore National Lab, + saransh1@llnl.gov >> @DATE: 10/28/2020 SS 1.0 original >> @DETAILS: constants for sphere sectors used for IPF coloring ''' @@ -263,7 +268,8 @@ def set_numba_cache(): 4. The NUMBA_CACHE_DIR environment variable is not defined If all of these are true, then numba will try to write to a - directory where it doesn't have permission, and cause the application to + directory where it doesn't have permission, + and cause the application to freeze. Avoid that by setting the cache dir ourselves. """ @@ -316,8 +322,10 @@ def is_writable_file(path): cRestmass = 9.1093837090E-31 ''' -adding another parametrization of the scattering factors. these are more recent -and more accurate. also used in Vesta (copied from there). see: +adding another parametrization of the +scattering factors. these are more recent +and more accurate. also used in Vesta +(copied from there). see: New Analytical coherent Scattering-Factor Functions for Free Atoms and Ions BY D. WAASMAIER AND A. KIRFEL @@ -638,7 +646,9 @@ def is_writable_file(path): 'Cf': ['0'] } ''' -this dictionary tabulates the small nuclear Thomson term fNT for all elements up to Z=92 +this dictionary tabulates the small +nuclear Thomson term fNT for all +elements up to Z=92 ''' fNT = { 'H':-0.00054423,'He':-0.00054817,'Li':-0.00071131,'Be':-0.00097394,'B':-0.0012687,'C':-0.0016442,'N':-0.0019191,'O':-0.0021944, @@ -737,8 +747,10 @@ def is_writable_file(path): 229 ]) -''' this variable encodes all the generators (including translations) for all 230 space groups - will be used to compute the full space group symmetry operators +''' this variable encodes all the generators +(including translations) for all 230 space groups +will be used to compute the full space group symmetry +operators ''' SYM_GL = [ "000 ", "100 ", "01cOOO0 ", @@ -996,8 +1008,10 @@ def is_writable_file(path): } ''' -this dictionary contains the generators encoded in each letter of the generator string -the full symmetry is generated by the repeated action of the generator matrix +this dictionary contains the generators encoded +in each letter of the generator string +the full symmetry is generated by the repeated +action of the generator matrix ''' ''' rotational, inversions, mirrors etc. components @@ -1100,49 +1114,21 @@ def is_writable_file(path): SYM_GENERATORS['Y'] = -1./4. SYM_GENERATORS['Z'] = -1./8. -''' this dictionary contains the space group numbers -which have two origin choices the two values are the -site symmetries of the origin. There are 24 such space -groups -''' -two_origin_choice = { -48 : ['222', '-1'], -50 : ['222/n', '-1'], -59 : ['mm2/n', '-1'], -68 : ['222', '-1'], -70 : ['222', '-1'], -85 : ['-4', '-1'], -86 : ['-4', '-1'], -88 : ['-4', '-1'], -125: ['422', '2/m'], -126: ['422/n', '-1'], -129: ['-4m2', '2/m'], -130: ['-4/ncn', '-1'], -133: ['-4121/c', '-1'], -134: ['-42m', '2/m'], -137: ['-4m2/n', '-1'], -138: ['-4cg', '2/m'], -141: ['-4m2', '2/m'], -142: ['-4c21', '-1'], -201: ['23', '-3'], -203: ['23', '-3'], -222: ['432', '-3'], -224: ['-43m', '-3m'], -227: ['-43m', '-3m'], -228: ['23', '-3'], -} - - ''' - @AUTHOR Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov + @AUTHOR Saransh Singh, + Lawrence Livermore National Lab, + saransh1@llnl.gov @DATE 11/23/2020 SS 1.0 original - @DETAIL. this list of symbols will help us to genrate the point group symmetries - in the cartesian space for any point group. this is needed for the - supergroup symmetry usd in the coloring scheme used in the package. this - needs to be a separate set of routines because the supergroup can be a - point group which is not the laue group of the crystal (e.g. m-3 --> m-3m) - the notation used will be the same as the one used for the space group - without any translations. + @DETAIL. this list of symbols will help us to genrate + the point group symmetries in the cartesian + space for any point group. this is needed for + the supergroup symmetry usd in the coloring + scheme used in the package. this needs to be a + separate set of routines because the supergroup + can be a point group which is not the laue group + of the crystal (e.g. m-3 --> m-3m) the notation + used will be the same as the one used for the + space group without any translations. ''' SYM_GL_PG = { 'c1': '1a', # only identity rotation diff --git a/hexrd/material/material.py b/hexrd/material/material.py index 5ca1ca67a..b447c6fd2 100644 --- a/hexrd/material/material.py +++ b/hexrd/material/material.py @@ -36,11 +36,11 @@ from hexrd.material.crystallography import PlaneData as PData from hexrd.material import symmetry, unitcell +from hexrd.material.symbols import two_origin_choice from hexrd.valunits import valWUnit from hexrd.constants import (ptable, - ptableinverse, - chargestate, - two_origin_choice) + ptableinverse, + chargestate) from os import path from pathlib import Path @@ -79,6 +79,7 @@ def _kev(x): def _key(x): return x.name + def get_default_sgsetting(sgnum): if sgnum in two_origin_choice: return 1 @@ -187,8 +188,8 @@ def __init__( elif isinstance(material_file, h5py.Group) or form in h5_suffixes: self._readHDFxtal(fhdf=material_file, xtal=name) if sgsetting is not None: - if sgsetting in [0, 1]: - self._sgsetting = sgsetting + if sgsetting in [0, 1]: + self._sgsetting = sgsetting else: # default name self._name = Material.DFLT_XTAL @@ -1436,7 +1437,7 @@ def loadMaterialList(cfgFile): def load_materials_hdf5( - f, dmin=None, kev=None): + f, dmin=None, kev=None): """Load materials from an HDF5 file The file uses the HDF5 file format. diff --git a/hexrd/material/mksupport.py b/hexrd/material/mksupport.py index 083272624..73f718700 100644 --- a/hexrd/material/mksupport.py +++ b/hexrd/material/mksupport.py @@ -1,6 +1,9 @@ -from hexrd.material.symbols import pstr_Elements, sitesym, \ - tworig, PrintPossibleSG, TRIG, pstr_spacegroup,\ - pstr_mkxtal +from hexrd.material.symbols import (pstr_Elements, + two_origin_choice, + PrintPossibleSG, + TRIG, + pstr_spacegroup, + pstr_mkxtal) import h5py import os import numpy as np @@ -8,6 +11,7 @@ import getpass from hexrd.material.unitcell import _StiffnessDict, _pgDict + def mk(filename, xtalname): # print some initial information for the user @@ -23,7 +27,7 @@ def mk(filename, xtalname): AtomInfo = GetAtomInfo() AtomInfo.update({'file': filename, 'xtalname': xtalname, - 'xtal_sys': xtal_sys, 'SG': space_group,\ + 'xtal_sys': xtal_sys, 'SG': space_group, 'SGsetting': iset}) Write2H5File(AtomInfo, lat_param) @@ -275,14 +279,15 @@ def GetSpaceGroup(xtal_sys, btrigonal, bhexset): def SpaceGroupSetting(sgnum): iset = 1 - if(sgnum in tworig): - idx = tworig.index(sgnum) + if(sgnum in two_origin_choice): + sitesym1 = two_origin_choice[sgnum][0] + sitesym2 = two_origin_choice[sgnum][1] print(' ---------------------------------------------') print(' This space group has two origin settings.') print(' The first setting has site symmetry : ' + - sitesym[2*idx - 2]) + sitesym1) print(' the second setting has site symmetry : ' + - sitesym[2*idx - 1]) + sitesym2) iset = input(' Which setting do you wish to use (1/2) : ') if(not iset.isdigit()): raise ValueError("Invalid integer value for atomic number.") @@ -294,7 +299,6 @@ def SpaceGroupSetting(sgnum): return iset-1 - def GetAtomInfo(): print(pstr_Elements) diff --git a/hexrd/material/symbols.py b/hexrd/material/symbols.py index ccc65ef36..18d8479b4 100644 --- a/hexrd/material/symbols.py +++ b/hexrd/material/symbols.py @@ -120,17 +120,37 @@ '89:Ac 90:Th 91:Pa 92:U' + "\n" \ ' ----------------------------------------------------------------------------------------------------------' -sitesym = ['222 ',' -1 ','222/n ',' -1 ','mm2/n ',' -1 ', \ -'222 ',' -1 ','222 ',' -1 ','-4 ',' -1 ', \ -'-4 ',' -1 ','-4 ',' -1 ','422 ','2/m ', \ -'422/n ',' -1 ','-4m2 ','2/m ','-4/ncn ',' -1 ', \ -'-4121/c',' -1 ','-42m ','2/m ','-4m2/n ',' -1 ', \ -'-4cg ','2/m ','-4m2 ','2/m ','-4c21 ',' -1 ', \ -'23 ',' -3 ','23 ',' -3 ','432 ',' -3 ', \ -'-43m ','-3m ','-43m ','-3m ','23 ',' -3 '] - -tworig = [48,50,59,68,70,85,86,88,125,126,129,130,133,134,137,138,\ -141,142,201,203,222,224,227,228] +''' this dictionary contains the space group numbers +which have two origin choices the two values are the +site symmetries of the origin. There are 24 such space +groups +''' +two_origin_choice = { +48 : ['222', '-1'], +50 : ['222/n', '-1'], +59 : ['mm2/n', '-1'], +68 : ['222', '-1'], +70 : ['222', '-1'], +85 : ['-4', '-1'], +86 : ['-4', '-1'], +88 : ['-4', '-1'], +125: ['422', '2/m'], +126: ['422/n', '-1'], +129: ['-4m2', '2/m'], +130: ['-4/ncn', '-1'], +133: ['-4121/c', '-1'], +134: ['-42m', '2/m'], +137: ['-4m2/n', '-1'], +138: ['-4cg', '2/m'], +141: ['-4m2', '2/m'], +142: ['-4c21', '-1'], +201: ['23', '-3'], +203: ['23', '-3'], +222: ['432', '-3'], +224: ['-43m', '-3m'], +227: ['-43m', '-3m'], +228: ['23', '-3'], +} def PrintPossibleSG(xtal_sys): diff --git a/hexrd/material/symmetry.py b/hexrd/material/symmetry.py index 5ea30a1bb..0d614f7ff 100644 --- a/hexrd/material/symmetry.py +++ b/hexrd/material/symmetry.py @@ -31,9 +31,9 @@ import numpy as np from numba import njit -from numpy import array, sqrt, pi, \ - vstack, c_, dot, \ - argmax +from numpy import (array, sqrt, pi, + vstack, c_, dot, + argmax) # from hexrd.rotations import quatOfAngleAxis, quatProductMatrix, fixQuat from hexrd import rotations as rot @@ -41,8 +41,9 @@ from hexrd.utils.decorators import memoize # Imports in case others are importing from here -from hexrd.rotations import toFundamentalRegion, ltypeOfLaueGroup, \ - quatOfLaueGroup +from hexrd.rotations import (toFundamentalRegion, + ltypeOfLaueGroup, + quatOfLaueGroup) # ============================================================================= @@ -78,8 +79,8 @@ def GeneratorString(sgnum): return constants.SYM_GL[sg] -def MakeGenerators(genstr, setting): +def MakeGenerators(genstr, setting): t = 'aOOO' mat = SYM_fillgen(t) @@ -99,7 +100,7 @@ def MakeGenerators(genstr, setting): if(n > 0): for i in range(n): istart = 2 + i * 4 - istop = 2 + (i+1) * 4 + istop = 2 + (i+1) * 4 t = genstr[istart:istop] @@ -114,39 +115,39 @@ def MakeGenerators(genstr, setting): if(genstr[istop] != '0'): if(setting != 0): t = genstr[istop+1:istop+4] - t = 'a' + t # get the translation without any rotation - sym = np.squeeze(SYM_fillgen(t, sgn=-1)) + t = 'a' + t # get the translation without any rotation + sym = np.squeeze(SYM_fillgen(t, sgn=-1)) sym2 = np.squeeze(SYM_fillgen(t)) for i in range(1, genmat.shape[0]): generator = np.dot(sym2, np.dot( - np.squeeze(genmat[i,:,:]), - sym)) - frac = np.modf(generator[0:3,3])[0] + np.squeeze(genmat[i, :, :]), + sym)) + frac = np.modf(generator[0:3, 3])[0] frac[frac < 0.] += 1. frac[np.abs(frac) < 1E-5] = 0.0 frac[np.abs(frac-1.0) < 1E-5] = 0.0 - generator[0:3,3] = frac - genmat[i,:,:] = generator + generator[0:3, 3] = frac + genmat[i, :, :] = generator return genmat, centrosymmetric + def SYM_fillgen(t, sgn=1): - mat = np.zeros([4,4]) - mat[3,3] = 1. + mat = np.zeros([4, 4]) + mat[3, 3] = 1. - mat[0:3,0:3] = constants.SYM_GENERATORS[t[0]] - mat[0:3,3] = sgn*np.array([constants.SYM_GENERATORS[t[1]], - constants.SYM_GENERATORS[t[2]], - constants.SYM_GENERATORS[t[3]] - ]) + mat[0:3, 0:3] = constants.SYM_GENERATORS[t[0]] + mat[0:3, 3] = sgn*np.array([constants.SYM_GENERATORS[t[1]], + constants.SYM_GENERATORS[t[2]], + constants.SYM_GENERATORS[t[3]] + ]) - mat = np.broadcast_to(mat, [1,4,4]) + mat = np.broadcast_to(mat, [1, 4, 4]) return mat @memoize(maxsize=20) def GenerateSGSym(sgnum, setting=0): - ''' get the generators for a space group using the generator string @@ -173,22 +174,22 @@ def GenerateSGSym(sgnum, setting=0): k1 = 0 while k1 < nsym: - g1 = np.squeeze(SYM_SG[k1,:,:]) + g1 = np.squeeze(SYM_SG[k1, :, :]) k2 = k1 while k2 < nsym: - g2 = np.squeeze(SYM_SG[k2,:,:]) + g2 = np.squeeze(SYM_SG[k2, :, :]) gnew = np.dot(g1, g2) # only fractional parts - frac = np.modf(gnew[0:3,3])[0] + frac = np.modf(gnew[0:3, 3])[0] frac[frac < 0.] += 1. frac[np.abs(frac) < 1E-5] = 0.0 frac[np.abs(frac-1.0) < 1E-5] = 0.0 - gnew[0:3,3] = frac + gnew[0:3, 3] = frac if(isnew(gnew, SYM_SG)): - gnew = np.broadcast_to(gnew, [1,4,4]) + gnew = np.broadcast_to(gnew, [1, 4, 4]) SYM_SG = np.concatenate((SYM_SG, gnew)) nsym += 1 @@ -203,12 +204,12 @@ def GenerateSGSym(sgnum, setting=0): SYM_PG_d_laue = GeneratePGSym_Laue(SYM_PG_d) for s in SYM_PG_d: - if(np.allclose(-np.eye(3),s)): + if(np.allclose(-np.eye(3), s)): centrosymmetric = True - return SYM_SG, SYM_PG_d, SYM_PG_d_laue, centrosymmetric, symmorphic + def GeneratePGSym(SYM_SG): ''' calculate the direct space point group symmetries @@ -222,19 +223,20 @@ def GeneratePGSym(SYM_SG): nsgsym = SYM_SG.shape[0] # first fill the identity rotation - SYM_PG_d = SYM_SG[0,0:3,0:3] - SYM_PG_d = np.broadcast_to(SYM_PG_d,[1,3,3]) - - for i in range(1,nsgsym): - g = SYM_SG[i,:,:] - t = g[0:3,3] - g = g[0:3,0:3] - if(isnew(g,SYM_PG_d)): - g = np.broadcast_to(g,[1,3,3]) + SYM_PG_d = SYM_SG[0, 0:3, 0:3] + SYM_PG_d = np.broadcast_to(SYM_PG_d, [1, 3, 3]) + + for i in range(1, nsgsym): + g = SYM_SG[i, :, :] + t = g[0:3, 3] + g = g[0:3, 0:3] + if(isnew(g, SYM_PG_d)): + g = np.broadcast_to(g, [1, 3, 3]) SYM_PG_d = np.concatenate((SYM_PG_d, g)) return SYM_PG_d.astype(np.int32) + def GeneratePGSym_Laue(SYM_PG_d): ''' generate the laue group symmetry for the given set of @@ -248,7 +250,7 @@ def GeneratePGSym_Laue(SYM_PG_d): first check if the group already has the inversion symmetry ''' for s in SYM_PG_d: - if(np.allclose(s,-np.eye(3))): + if(np.allclose(s, -np.eye(3))): return SYM_PG_d ''' @@ -256,7 +258,7 @@ def GeneratePGSym_Laue(SYM_PG_d): add the inversion symmetry ''' SYM_PG_d_laue = SYM_PG_d - g = np.broadcast_to(-np.eye(3).astype(np.int32),[1,3,3]) + g = np.broadcast_to(-np.eye(3).astype(np.int32), [1, 3, 3]) SYM_PG_d_laue = np.concatenate((SYM_PG_d_laue, g)) ''' @@ -266,14 +268,14 @@ def GeneratePGSym_Laue(SYM_PG_d): nsym = SYM_PG_d_laue.shape[0] k1 = 0 while k1 < nsym: - g1 = np.squeeze(SYM_PG_d_laue[k1,:,:]) + g1 = np.squeeze(SYM_PG_d_laue[k1, :, :]) k2 = k1 while k2 < nsym: - g2 = np.squeeze(SYM_PG_d_laue[k2,:,:]) + g2 = np.squeeze(SYM_PG_d_laue[k2, :, :]) gnew = np.dot(g1, g2) if(isnew(gnew, SYM_PG_d_laue)): - gnew = np.broadcast_to(gnew, [1,3,3]) + gnew = np.broadcast_to(gnew, [1, 3, 3]) SYM_PG_d_laue = np.concatenate((SYM_PG_d_laue, gnew)) nsym += 1 @@ -295,6 +297,7 @@ def isnew(mat, sym_mats): return False return True + def latticeType(sgnum): if(sgnum <= 2): @@ -309,7 +312,7 @@ def latticeType(sgnum): return 'trigonal' elif(sgnum > 167 and sgnum <= 194): return 'hexagonal' - elif(sgnum > 194 and sgnum <=230): + elif(sgnum > 194 and sgnum <= 230): return 'cubic' else: raise RuntimeError('symmetry.latticeType: unknown space group number') @@ -323,14 +326,15 @@ def MakeGenerators_PGSYM(pggenstr): for any point group. this is needed for the coloring routines ''' ngen = int(pggenstr[0]) - SYM_GEN_PG = np.zeros([ngen,3,3]) + SYM_GEN_PG = np.zeros([ngen, 3, 3]) for i in range(ngen): s = pggenstr[i+1] - SYM_GEN_PG[i,:,:] = constants.SYM_GENERATORS[s] + SYM_GEN_PG[i, :, :] = constants.SYM_GENERATORS[s] return SYM_GEN_PG + def GeneratePGSYM(pgsym): ''' @AUTHOR Saransh Singh, Lawrence Livermore National Lab, saransh1@llnl.gov @@ -345,7 +349,6 @@ def GeneratePGSYM(pgsym): generate the powers of the group ''' - ''' now go through the group actions and see if its a new matrix if it is then add it to the group @@ -353,14 +356,14 @@ def GeneratePGSYM(pgsym): nsym = SYM_GEN_PG.shape[0] k1 = 0 while k1 < nsym: - g1 = np.squeeze(SYM_GEN_PG[k1,:,:]) + g1 = np.squeeze(SYM_GEN_PG[k1, :, :]) k2 = k1 while k2 < nsym: - g2 = np.squeeze(SYM_GEN_PG[k2,:,:]) + g2 = np.squeeze(SYM_GEN_PG[k2, :, :]) gnew = np.dot(g1, g2) if(isnew(gnew, SYM_GEN_PG)): - gnew = np.broadcast_to(gnew, [1,3,3]) + gnew = np.broadcast_to(gnew, [1, 3, 3]) SYM_GEN_PG = np.concatenate((SYM_GEN_PG, gnew)) nsym += 1 From f0f6c1b3da1ff7b08ab45390bd5f23a10df13e41 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Tue, 20 Aug 2024 11:46:44 -0700 Subject: [PATCH 11/16] more PEP8 --- hexrd/material/symbols.py | 278 ++++++++++++++++++++------------------ 1 file changed, 144 insertions(+), 134 deletions(-) diff --git a/hexrd/material/symbols.py b/hexrd/material/symbols.py index 18d8479b4..200ae37ca 100644 --- a/hexrd/material/symbols.py +++ b/hexrd/material/symbols.py @@ -1,6 +1,5 @@ pstr_mkxtal = "\n\n This is a program to create a HDF5 file for storing crystallographic information.\n " -pstr_mkxtal = pstr_mkxtal + " This format is the same format as used in the EMsoft (electron microscoy) suite.\n " pstr_mkxtal = pstr_mkxtal + " The following inputs are required:\n " pstr_mkxtal = pstr_mkxtal + " Crystal System:\n" pstr_mkxtal = pstr_mkxtal + " 1. Cubic\n" @@ -24,88 +23,98 @@ pstr_mkxtal = pstr_mkxtal + " trigonal symmetry is the hexagonal setting. When you select\n" pstr_mkxtal = pstr_mkxtal + " crystal system 5 above, you will be prompted for the setting. \n" -pstr_spacegroup = [ " P 1 " ," P -1 ", \ -# MONOCLINIC SPACE GROUPS -" P 2 " ," P 21 " ," C 2 " ," P m ", \ -" P c " ," C m " ," C c " ," P 2/m ", \ -" P 21/m " ," C 2/m " ," P 2/c " ," P 21/c ", \ -" C 2/c ", \ -# ORTHORHOMBIC SPACE GROUPS -" P 2 2 2 " ," P 2 2 21 " ," P 21 21 2 " ," P 21 21 21", \ -" C 2 2 21 " ," C 2 2 2 " ," F 2 2 2 " ," I 2 2 2 ", \ -" I 21 21 21" ," P m m 2 " ," P m c 21 " ," P c c 2 ", \ -" P m a 2 " ," P c a 21 " ," P n c 2 " ," P m n 21 ", \ -" P b a 2 " ," P n a 21 " ," P n n 2 " ," C m m 2 ", \ -" C m c 21 " ," C c c 2 " ," A m m 2 " ," A b m 2 ", \ -" A m a 2 " ," A b a 2 " ," F m m 2 " ," F d d 2 ", \ -" I m m 2 " ," I b a 2 " ," I m a 2 " ," P m m m ", \ -" P n n n " ," P c c m " ," P b a n " ," P m m a ", \ -" P n n a " ," P m n a " ," P c c a " ," P b a m ", \ -" P c c n " ," P b c m " ," P n n m " ," P m m n ", \ -" P b c n " ," P b c a " ," P n m a " ," C m c m ", \ -" C m c a " ," C m m m " ," C c c m " ," C m m a ", \ -" C c c a " ," F m m m " ," F d d d " ," I m m m ", \ -" I b a m " ," I b c a " ," I m m a ", \ -# TETRAGONAL SPACE GROUPS -" P 4 " ," P 41 " ," P 42 " ," P 43 ", \ -" I 4 " ," I 41 " ," P -4 " ," I -4 ", \ -" P 4/m " ," P 42/m " ," P 4/n " ," P 42/n ", \ -" I 4/m " ," I 41/a " ," P 4 2 2 " ," P 4 21 2 ", \ -" P 41 2 2 " ," P 41 21 2 " ," P 42 2 2 " ," P 42 21 2 ", \ -" P 43 2 2 " ," P 43 21 2 " ," I 4 2 2 " ," I 41 2 2 ", \ -" P 4 m m " ," P 4 b m " ," P 42 c m " ," P 42 n m ", \ -" P 4 c c " ," P 4 n c " ," P 42 m c " ," P 42 b c ", \ -" I 4 m m " ," I 4 c m " ," I 41 m d " ," I 41 c d ", \ -" P -4 2 m " ," P -4 2 c " ," P -4 21 m " ," P -4 21 c ", \ -" P -4 m 2 " ," P -4 c 2 " ," P -4 b 2 " ," P -4 n 2 ", \ -" I -4 m 2 " ," I -4 c 2 " ," I -4 2 m " ," I -4 2 d ", \ -" P 4/m m m " ," P 4/m c c " ," P 4/n b m " ," P 4/n n c ", \ -" P 4/m b m " ," P 4/m n c " ," P 4/n m m " ," P 4/n c c ", \ -" P 42/m m c" ," P 42/m c m" ," P 42/n b c" ," P 42/n n m", \ -" P 42/m b c" ," P 42/m n m" ," P 42/n m c" ," P 42/n c m", \ -" I 4/m m m " ," I 4/m c m " ," I 41/a m d" ," I 41/a c d", \ -# RHOMBOHEDRAL SPACE GROUPS -" P 3 " ," P 31 " ," P 32 " ," R 3 ", \ -" P -3 " ," R -3 " ," P 3 1 2 " ," P 3 2 1 ", \ -" P 31 1 2 " ," P 31 2 1 " ," P 32 1 2 " ," P 32 2 1 ", \ -" R 3 2 " ," P 3 m 1 " ," P 3 1 m " ," P 3 c 1 ", \ -" P 3 1 c " ," R 3 m " ," R 3 c " ," P -3 1 m ", \ -" P -3 1 c " ," P -3 m 1 " ," P -3 c 1 " ," R -3 m ", \ -" R -3 c ", \ -# HEXAGONAL SPACE GROUPS -" P 6 " ," P 61 " ," P 65 " ," P 62 ", \ -" P 64 " ," P 63 " ," P -6 " ," P 6/m ", \ -" P 63/m " ," P 6 2 2 " ," P 61 2 2 " ," P 65 2 2 ", \ -" P 62 2 2 " ," P 64 2 2 " ," P 63 2 2 " ," P 6 m m ", \ -" P 6 c c " ," P 63 c m " ," P 63 m c " ," P -6 m 2 ", \ -" P -6 c 2 " ," P -6 2 m " ," P -6 2 c " ," P 6/m m m ", \ -" P 6/m c c " ," P 63/m c m" ," P 63/m m c", \ -#CUBIC SPACE GROUPS -" P 2 3 " ," F 2 3 " ," I 2 3 " ," P 21 3 ", \ -" I 21 3 " ," P m 3 " ," P n 3 " ," F m 3 ", \ -" F d 3 " ," I m 3 " ," P a 3 " ," I a 3 ", \ -" P 4 3 2 " ," P 42 3 2 " ," F 4 3 2 " ," F 41 3 2 ", \ -" I 4 3 2 " ," P 43 3 2 " ," P 41 3 2 " ," I 41 3 2 ", \ -" P -4 3 m " ," F -4 3 m " ," I -4 3 m " ," P -4 3 n ", \ -" F -4 3 c " ," I -4 3 d " ," P m 3 m " ," P n 3 n ", \ -" P m 3 n " ," P n 3 m " ," F m 3 m " ," F m 3 c ", \ -" F d 3 m " ," F d 3 c " ," I m 3 m " ," I a 3 d ", \ -# TRIGONAL GROUPS RHOMBOHEDRAL SETTING -" R 3 |146" ," R -3 |148" ," R 3 2 |155" ," R 3 m |160", \ -" R 3 c |161" ," R -3 m|166" ," R -3 c|167"] +pstr_spacegroup = +[ + " P 1 ", " P -1 ", \ + # MONOCLINIC SPACE GROUPS + " P 2 ", " P 21 ", " C 2 ", " P m ", \ + " P c ", " C m ", " C c ", " P 2/m ", \ + " P 21/m ", " C 2/m ", " P 2/c ", " P 21/c ", \ + " C 2/c ", \ + # ORTHORHOMBIC SPACE GROUPS + " P 2 2 2 ", " P 2 2 21 ", " P 21 21 2 ", " P 21 21 21", \ + " C 2 2 21 ", " C 2 2 2 ", " F 2 2 2 ", " I 2 2 2 ", \ + " I 21 21 21", " P m m 2 ", " P m c 21 ", " P c c 2 ", \ + " P m a 2 ", " P c a 21 ", " P n c 2 ", " P m n 21 ", \ + " P b a 2 ", " P n a 21 ", " P n n 2 ", " C m m 2 ", \ + " C m c 21 ", " C c c 2 ", " A m m 2 ", " A b m 2 ", \ + " A m a 2 ", " A b a 2 ", " F m m 2 ", " F d d 2 ", \ + " I m m 2 ", " I b a 2 ", " I m a 2 ", " P m m m ", \ + " P n n n ", " P c c m ", " P b a n ", " P m m a ", \ + " P n n a ", " P m n a ", " P c c a ", " P b a m ", \ + " P c c n ", " P b c m ", " P n n m ", " P m m n ", \ + " P b c n ", " P b c a ", " P n m a ", " C m c m ", \ + " C m c a ", " C m m m ", " C c c m ", " C m m a ", \ + " C c c a ", " F m m m ", " F d d d ", " I m m m ", \ + " I b a m ", " I b c a ", " I m m a ", \ + # TETRAGONAL SPACE GROUPS + " P 4 ", " P 41 ", " P 42 ", " P 43 ", \ + " I 4 ", " I 41 ", " P -4 ", " I -4 ", \ + " P 4/m ", " P 42/m ", " P 4/n ", " P 42/n ", \ + " I 4/m ", " I 41/a ", " P 4 2 2 ", " P 4 21 2 ", \ + " P 41 2 2 ", " P 41 21 2 ", " P 42 2 2 ", " P 42 21 2 ", \ + " P 43 2 2 ", " P 43 21 2 ", " I 4 2 2 ", " I 41 2 2 ", \ + " P 4 m m ", " P 4 b m ", " P 42 c m ", " P 42 n m ", \ + " P 4 c c ", " P 4 n c ", " P 42 m c ", " P 42 b c ", \ + " I 4 m m ", " I 4 c m ", " I 41 m d ", " I 41 c d ", \ + " P -4 2 m ", " P -4 2 c ", " P -4 21 m ", " P -4 21 c ", \ + " P -4 m 2 ", " P -4 c 2 ", " P -4 b 2 ", " P -4 n 2 ", \ + " I -4 m 2 ", " I -4 c 2 ", " I -4 2 m ", " I -4 2 d ", \ + " P 4/m m m ", " P 4/m c c ", " P 4/n b m ", " P 4/n n c ", \ + " P 4/m b m ", " P 4/m n c ", " P 4/n m m ", " P 4/n c c ", \ + " P 42/m m c", " P 42/m c m", " P 42/n b c", " P 42/n n m", \ + " P 42/m b c", " P 42/m n m", " P 42/n m c", " P 42/n c m", \ + " I 4/m m m ", " I 4/m c m ", " I 41/a m d", " I 41/a c d", \ + # RHOMBOHEDRAL SPACE GROUPS + " P 3 ", " P 31 ", " P 32 ", " R 3 ", \ + " P -3 ", " R -3 ", " P 3 1 2 ", " P 3 2 1 ", \ + " P 31 1 2 ", " P 31 2 1 ", " P 32 1 2 ", " P 32 2 1 ", \ + " R 3 2 ", " P 3 m 1 ", " P 3 1 m ", " P 3 c 1 ", \ + " P 3 1 c ", " R 3 m ", " R 3 c ", " P -3 1 m ", \ + " P -3 1 c ", " P -3 m 1 ", " P -3 c 1 ", " R -3 m ", \ + " R -3 c ", \ + # HEXAGONAL SPACE GROUPS + " P 6 ", " P 61 ", " P 65 ", " P 62 ", \ + " P 64 ", " P 63 ", " P -6 ", " P 6/m ", \ + " P 63/m ", " P 6 2 2 ", " P 61 2 2 ", " P 65 2 2 ", \ + " P 62 2 2 ", " P 64 2 2 ", " P 63 2 2 ", " P 6 m m ", \ + " P 6 c c ", " P 63 c m ", " P 63 m c ", " P -6 m 2 ", \ + " P -6 c 2 ", " P -6 2 m ", " P -6 2 c ", " P 6/m m m ", \ + " P 6/m c c ", " P 63/m c m", " P 63/m m c", \ + # CUBIC SPACE GROUPS + " P 2 3 ", " F 2 3 ", " I 2 3 ", " P 21 3 ", \ + " I 21 3 ", " P m 3 ", " P n 3 ", " F m 3 ", \ + " F d 3 ", " I m 3 ", " P a 3 ", " I a 3 ", \ + " P 4 3 2 ", " P 42 3 2 ", " F 4 3 2 ", " F 41 3 2 ", \ + " I 4 3 2 ", " P 43 3 2 ", " P 41 3 2 ", " I 41 3 2 ", \ + " P -4 3 m ", " F -4 3 m ", " I -4 3 m ", " P -4 3 n ", \ + " F -4 3 c ", " I -4 3 d ", " P m 3 m ", " P n 3 n ", \ + " P m 3 n ", " P n 3 m ", " F m 3 m ", " F m 3 c ", \ + " F d 3 m ", " F d 3 c ", " I m 3 m ", " I a 3 d ", \ + # TRIGONAL GROUPS RHOMBOHEDRAL SETTING + " R 3 |146", " R -3 |148", " R 3 2 |155", " R 3 m |160", \ + " R 3 c |161", " R -3 m|166", " R -3 c|167" +] -xtal_dict = {1:'cubic', 2:'tetragonal', 3:'orthorhombic', 4:'hexagonal', 5:'trigonal', 6:'monoclinic', 7:'triclinic'} -xtal_sys_dict = {'cubic':1, 'tetragonal':2, 'orthorhombic':3, 'hexagonal':4, 'trigonal':5, 'monoclinic':6, 'triclinic':7} +xtal_dict = {1: 'cubic', 2: 'tetragonal', 3: 'orthorhombic', + 4: 'hexagonal', 5: 'trigonal', 6: 'monoclinic', + 7: 'triclinic'} +xtal_sys_dict = {'cubic': 1, 'tetragonal': 2, 'orthorhombic': 3, + 'hexagonal': 4, 'trigonal': 5, 'monoclinic': 6, + 'triclinic': 7} -pstr_pointgroup = [' 1',' -1',' 2',' m',' 2/m',' 222', \ -' mm2',' mmm',' 4',' -4',' 4/m',' 422', \ -' 4mm',' -42m','4/mmm',' 3',' -3',' 32', \ -' 3m',' -3m',' 6',' -6',' 6/m',' 622', \ -' 6mm',' -6m2','6/mmm',' 23',' m3',' 432', \ -' -43m',' m-3m',' 532',' 822',' 1022',' 1222' ] -TRIG = [146,148,155,160,161,166,167] + +pstr_pointgroup = [' 1', ' -1', ' 2', ' m', ' 2/m', ' 222', + ' mm2', ' mmm', ' 4', ' -4', ' 4/m', ' 422', + ' 4mm', ' -42m', '4/mmm', ' 3', ' -3', ' 32', + ' 3m', ' -3m', ' 6', ' -6', ' 6/m', ' 622', + ' 6mm', ' -6m2', '6/mmm', ' 23', ' m3', ' 432', + ' -43m', ' m-3m', ' 532', ' 822', ' 1022', ' 1222'] + +TRIG = [146, 148, 155, 160, 161, 166, 167] + # symbols and Z for all elements pstr_Elements = ' ------------------------------------ Periodic Table of the Elements --------------------------------------' + "\n" \ @@ -126,66 +135,66 @@ groups ''' two_origin_choice = { -48 : ['222', '-1'], -50 : ['222/n', '-1'], -59 : ['mm2/n', '-1'], -68 : ['222', '-1'], -70 : ['222', '-1'], -85 : ['-4', '-1'], -86 : ['-4', '-1'], -88 : ['-4', '-1'], -125: ['422', '2/m'], -126: ['422/n', '-1'], -129: ['-4m2', '2/m'], -130: ['-4/ncn', '-1'], -133: ['-4121/c', '-1'], -134: ['-42m', '2/m'], -137: ['-4m2/n', '-1'], -138: ['-4cg', '2/m'], -141: ['-4m2', '2/m'], -142: ['-4c21', '-1'], -201: ['23', '-3'], -203: ['23', '-3'], -222: ['432', '-3'], -224: ['-43m', '-3m'], -227: ['-43m', '-3m'], -228: ['23', '-3'], + 48: ['222', '-1'], + 50: ['222/n', '-1'], + 59: ['mm2/n', '-1'], + 68: ['222', '-1'], + 70: ['222', '-1'], + 85: ['-4', '-1'], + 86: ['-4', '-1'], + 88: ['-4', '-1'], + 125: ['422', '2/m'], + 126: ['422/n', '-1'], + 129: ['-4m2', '2/m'], + 130: ['-4/ncn', '-1'], + 133: ['-4121/c', '-1'], + 134: ['-42m', '2/m'], + 137: ['-4m2/n', '-1'], + 138: ['-4cg', '2/m'], + 141: ['-4m2', '2/m'], + 142: ['-4c21', '-1'], + 201: ['23', '-3'], + 203: ['23', '-3'], + 222: ['432', '-3'], + 224: ['-43m', '-3m'], + 227: ['-43m', '-3m'], + 228: ['23', '-3'], } def PrintPossibleSG(xtal_sys): - if(xtal_sys == 1): - sgmax = 230 - sgmin = 195 - elif(xtal_sys == 2): - sgmax = 142 - sgmin = 75 - elif(xtal_sys == 3): - sgmax = 74 - sgmin = 16 - elif(xtal_sys == 4): - sgmax = 194 - sgmin = 168 - elif(xtal_sys == 5): - sgmax = 167 - sgmin = 143 - elif(xtal_sys == 6): - sgmax = 15 - sgmin = 3 - elif(xtal_sys == 7): - sgmax = 2 - sgmin = 1 + if(xtal_sys == 1): + sgmax = 230 + sgmin = 195 + elif(xtal_sys == 2): + sgmax = 142 + sgmin = 75 + elif(xtal_sys == 3): + sgmax = 74 + sgmin = 16 + elif(xtal_sys == 4): + sgmax = 194 + sgmin = 168 + elif(xtal_sys == 5): + sgmax = 167 + sgmin = 143 + elif(xtal_sys == 6): + sgmax = 15 + sgmin = 3 + elif(xtal_sys == 7): + sgmax = 2 + sgmin = 1 - for i in range(sgmin,sgmax+1): - j = i - sgmin + 1 - pstr = str(i) + ":" + pstr_spacegroup[i-1] + "\t" - if(j % 4 == 0 or j == sgmax): - print(pstr) - else: - print(pstr, end='') - print("\n") + for i in range(sgmin, sgmax+1): + j = i - sgmin + 1 + pstr = str(i) + ":" + pstr_spacegroup[i-1] + "\t" + if(j % 4 == 0 or j == sgmax): + print(pstr) + else: + print(pstr, end='') + print("\n") - return sgmin, sgmax + return sgmin, sgmax # @@ -1257,6 +1266,7 @@ def PrintPossibleSG(xtal_sys): 230 I a -3 d """ + def _buildDict(hstr): """build the dictionaries from the notation string @@ -1267,7 +1277,7 @@ def _buildDict(hstr): for the rhombohedral lattices so that they use the hexagonal convention. """ - d = dict() + d = dict() di = dict() hs = hstr.split('\n') for l in hs: From 49432891c374aa85e8dbcd9d6269a151f1bcfdd2 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Tue, 20 Aug 2024 11:48:16 -0700 Subject: [PATCH 12/16] final few PEP8 --- hexrd/material/symbols.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/hexrd/material/symbols.py b/hexrd/material/symbols.py index 200ae37ca..a2a2cae8a 100644 --- a/hexrd/material/symbols.py +++ b/hexrd/material/symbols.py @@ -106,12 +106,16 @@ -pstr_pointgroup = [' 1', ' -1', ' 2', ' m', ' 2/m', ' 222', - ' mm2', ' mmm', ' 4', ' -4', ' 4/m', ' 422', - ' 4mm', ' -42m', '4/mmm', ' 3', ' -3', ' 32', - ' 3m', ' -3m', ' 6', ' -6', ' 6/m', ' 622', - ' 6mm', ' -6m2', '6/mmm', ' 23', ' m3', ' 432', - ' -43m', ' m-3m', ' 532', ' 822', ' 1022', ' 1222'] +pstr_pointgroup = +[ + ' 1', ' -1', ' 2', ' m', ' 2/m', ' 222', + ' mm2', ' mmm', ' 4', ' -4', ' 4/m', ' 422', + ' 4mm', ' -42m', '4/mmm', ' 3', ' -3', ' 32', + ' 3m', ' -3m', ' 6', ' -6', ' 6/m', ' 622', + ' 6mm', ' -6m2', '6/mmm', ' 23', ' m3', ' 432', + ' -43m', ' m-3m', ' 532', ' 822', ' 1022', ' 1222' +] + TRIG = [146, 148, 155, 160, 161, 166, 167] @@ -129,11 +133,12 @@ '89:Ac 90:Th 91:Pa 92:U' + "\n" \ ' ----------------------------------------------------------------------------------------------------------' -''' this dictionary contains the space group numbers +'''this dictionary contains the space group numbers which have two origin choices the two values are the site symmetries of the origin. There are 24 such space -groups -''' +groups''' + + two_origin_choice = { 48: ['222', '-1'], 50: ['222/n', '-1'], From 0d86fdc8b0856088c659a0bb75f82e061876a340 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Tue, 20 Aug 2024 11:52:01 -0700 Subject: [PATCH 13/16] accidentally introduced a typo --- hexrd/material/symbols.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/hexrd/material/symbols.py b/hexrd/material/symbols.py index a2a2cae8a..637a5136a 100644 --- a/hexrd/material/symbols.py +++ b/hexrd/material/symbols.py @@ -23,8 +23,7 @@ pstr_mkxtal = pstr_mkxtal + " trigonal symmetry is the hexagonal setting. When you select\n" pstr_mkxtal = pstr_mkxtal + " crystal system 5 above, you will be prompted for the setting. \n" -pstr_spacegroup = -[ +pstr_spacegroup = [ " P 1 ", " P -1 ", \ # MONOCLINIC SPACE GROUPS " P 2 ", " P 21 ", " C 2 ", " P m ", \ @@ -93,8 +92,8 @@ " F d 3 m ", " F d 3 c ", " I m 3 m ", " I a 3 d ", \ # TRIGONAL GROUPS RHOMBOHEDRAL SETTING " R 3 |146", " R -3 |148", " R 3 2 |155", " R 3 m |160", \ - " R 3 c |161", " R -3 m|166", " R -3 c|167" -] + " R 3 c |161", " R -3 m|166", " R -3 c|167"] + xtal_dict = {1: 'cubic', 2: 'tetragonal', 3: 'orthorhombic', @@ -106,15 +105,14 @@ -pstr_pointgroup = -[ +pstr_pointgroup = [ ' 1', ' -1', ' 2', ' m', ' 2/m', ' 222', ' mm2', ' mmm', ' 4', ' -4', ' 4/m', ' 422', ' 4mm', ' -42m', '4/mmm', ' 3', ' -3', ' 32', ' 3m', ' -3m', ' 6', ' -6', ' 6/m', ' 622', ' 6mm', ' -6m2', '6/mmm', ' 23', ' m3', ' 432', - ' -43m', ' m-3m', ' 532', ' 822', ' 1022', ' 1222' -] + ' -43m', ' m-3m', ' 532', ' 822', ' 1022', ' 1222'] + TRIG = [146, 148, 155, 160, 161, 166, 167] From ae5bb6b2eba1ddfbf5e72270053ed1a964f84b6c Mon Sep 17 00:00:00 2001 From: Zack Date: Tue, 20 Aug 2024 15:05:38 -0400 Subject: [PATCH 14/16] Resolve compiler flags for Windows (#706) * Resolve compiler flags for Windows * Remove cpp std identifier * Update pyproject * Re-add std c++14 * Specify c++14 only for c++, not c * Remove erroneous c++ std annotation * Remove changes to setuptools_scm version scheme These were not actually necessary to fix the mac packaging after all. Signed-off-by: Patrick Avery --------- Signed-off-by: Patrick Avery Co-authored-by: Patrick Avery --- setup.py | 48 ++++++++++++++++++------------------------------ 1 file changed, 18 insertions(+), 30 deletions(-) diff --git a/setup.py b/setup.py index 88a5b63e1..ae91a9620 100644 --- a/setup.py +++ b/setup.py @@ -34,37 +34,36 @@ # Add it to the requirements automatically for intel computers. install_reqs.append('tbb') -# This will determine which compiler is being used to build the C modules +# Determine which compiler is being used to build the C/C++ modules compiler_type = distutils.ccompiler.get_default_compiler() if compiler_type in ("unix", "mingw32"): - compiler_optimize_flags = ['-O3', '-ftree-vectorize'] + compiler_flags = ['-O3', '-ftree-vectorize', '-Wall', '-funroll-loops'] + if not sys.platform.startswith('win'): + compiler_flags.append('-fPIC') elif compiler_type == "msvc": - compiler_optimize_flags = ['/Ox', '/GL'] + compiler_flags = ['/Ox', '/GL', '/std:c++14'] else: - compiler_optimize_flags = [] - + compiler_flags = [] -# extension for convolution from astropy +# Extension for convolution from astropy def get_convolution_extensions(): c_convolve_pkgdir = Path('hexrd') / 'convolution' src_files = [str(c_convolve_pkgdir / 'src/convolve.c')] - extra_compile_args = ['-UNDEBUG'] - if not sys.platform.startswith('win'): - extra_compile_args.append('-fPIC') - extra_compile_args += compiler_optimize_flags + extra_compile_args = ['-UNDEBUG'] + compiler_flags # Add '-Rpass-missed=.*' to ``extra_compile_args`` when compiling with # clang to report missed optimizations - _convolve_ext = Extension(name='hexrd.convolution._convolve', - sources=src_files, - extra_compile_args=extra_compile_args, - include_dirs=[numpy.get_include()], - language='c') + _convolve_ext = Extension( + name='hexrd.convolution._convolve', + sources=src_files, + extra_compile_args=extra_compile_args, + include_dirs=[numpy.get_include()], + language='c' + ) return [_convolve_ext] - def get_include_path(library_name): env_var_hint = os.getenv(f"{library_name.upper()}_INCLUDE_DIR") if env_var_hint is not None and os.path.exists(env_var_hint): @@ -97,7 +96,6 @@ def get_include_path(library_name): # It should exist now return full_path - def get_pybind11_include_path(): # If we can import pybind11, use that include path try: @@ -110,16 +108,10 @@ def get_pybind11_include_path(): # Otherwise, we will download the source and include that return get_include_path('pybind11') - def get_cpp_extensions(): cpp_transform_pkgdir = Path('hexrd') / 'transforms/cpp_sublibrary' src_files = [str(cpp_transform_pkgdir / 'src/inverse_distortion.cpp')] - extra_compile_args = ['-O3', '-Wall', '-shared', '-std=c++14', - '-funroll-loops'] - if not sys.platform.startswith('win'): - extra_compile_args.append('-fPIC') - # Define include directories include_dirs = [ get_include_path('xsimd'), @@ -131,14 +123,13 @@ def get_cpp_extensions(): inverse_distortion_ext = Extension( name='hexrd.extensions.inverse_distortion', sources=src_files, - extra_compile_args=extra_compile_args, + extra_compile_args=compiler_flags+['-std=c++14'], include_dirs=include_dirs, language='c++', ) return [inverse_distortion_ext] - def get_old_xfcapi_extension_modules(): # for transforms srclist = ['transforms_CAPI.c', 'transforms_CFUNC.c'] @@ -147,23 +138,21 @@ def get_old_xfcapi_extension_modules(): 'hexrd.extensions._transforms_CAPI', sources=srclist, include_dirs=[np_include_dir], - extra_compile_args=compiler_optimize_flags, + extra_compile_args=compiler_flags, ) return [transforms_mod] - def get_new_xfcapi_extension_modules(): transforms_mod = Extension( 'hexrd.extensions._new_transforms_capi', sources=['hexrd/transforms/new_capi/module.c'], include_dirs=[np_include_dir], - extra_compile_args=compiler_optimize_flags, + extra_compile_args=compiler_flags, ) return [transforms_mod] - def get_extension_modules(): # Flatten the lists return [item for sublist in ( @@ -173,7 +162,6 @@ def get_extension_modules(): get_cpp_extensions(), ) for item in sublist] - ext_modules = get_extension_modules() # use entry_points, not scripts: From 79c89837b306e5413666ab13cf1e3c67b1e1c9e2 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Tue, 20 Aug 2024 13:06:21 -0700 Subject: [PATCH 15/16] only select setting 2 by default if crystal is cubic e.g. diamond, silicon etc. --- hexrd/material/material.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/hexrd/material/material.py b/hexrd/material/material.py index b447c6fd2..df85eee6a 100644 --- a/hexrd/material/material.py +++ b/hexrd/material/material.py @@ -81,10 +81,13 @@ def _key(x): def get_default_sgsetting(sgnum): + setting = 0 if sgnum in two_origin_choice: - return 1 - else: - return 0 + # only set cubic materials + # in the second setting by default + if sgnum > 194: + setting = 1 + return setting # # ---------------------------------------------------CLASS: Material From 3601b9f32a6bbaa65caeaad6e229b5f612964108 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Tue, 20 Aug 2024 13:43:14 -0700 Subject: [PATCH 16/16] roll back changes to default space group setting, and modify tests to only assume setting 1 --- hexrd/material/material.py | 9 +++------ tests/planedata/test_with_data.py | 3 ++- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/hexrd/material/material.py b/hexrd/material/material.py index df85eee6a..b447c6fd2 100644 --- a/hexrd/material/material.py +++ b/hexrd/material/material.py @@ -81,13 +81,10 @@ def _key(x): def get_default_sgsetting(sgnum): - setting = 0 if sgnum in two_origin_choice: - # only set cubic materials - # in the second setting by default - if sgnum > 194: - setting = 1 - return setting + return 1 + else: + return 0 # # ---------------------------------------------------CLASS: Material diff --git a/tests/planedata/test_with_data.py b/tests/planedata/test_with_data.py index b4e771792..0b4e9d539 100644 --- a/tests/planedata/test_with_data.py +++ b/tests/planedata/test_with_data.py @@ -24,7 +24,8 @@ def materials(test_data_dir): for mat_name in material_names: # Load {test_data_dir}/materials/{mat_name}.cif mat = Material( - mat_name, str(test_data_dir) + "/materials/" + mat_name + ".cif" + mat_name, str(test_data_dir) + "/materials/" + mat_name + ".cif", + sgsetting=0 ) mats[mat_name] = mat.planeData return mats