Merge remote-tracking branch 'origin/master'

# Conflicts:
#	setup.cfg

.gitignore

@@ -12,6 +12,9 @@ __pycache__/
 # Protein DataBank cache folder
 pdb_cache/
 
+# Jupyter notebooks
+notebooks/
+
 # Distribution / packaging
 .Python
 env/

setup.py

@@ -6,6 +6,15 @@
 from setuptools import setup, find_packages
 #from Cython.Build import cythonize
 
+# How to submit a package on PyPi:
+# http://peterdowns.com/posts/first-time-with-pypi.html
+#
+# Test upload:
+#   > python setup.py sdist upload -r pypitest
+#
+# Real upload:
+#   > python setup.py sdist upload -r pypi
+
 BASE_PACKAGE = 'skued'
 
 wavelets = chain.from_iterable([glob('skued\\baseline\\data\\*.npy'), 

skued/__init__.py

@@ -4,11 +4,12 @@
 __license__ = 'MIT'
 __version__ = '0.4.4.1' # TODO: automatic versioning?
 
-from .array_utils import repeated_array
+from .affine import (affine_map, change_basis_mesh, change_of_basis, is_basis,
+                     is_rotation_matrix, minimum_image_distance,
+                     rotation_matrix, transform, translation_matrix,
+                     translation_rotation_matrix)
+from .array_utils import mirror, repeated_array
 from .parallel import pmap, preduce
 from .plot_utils import spectrum_colors
-from .quantities import lorentz, electron_wavelength, interaction_parameter
-from .affine import (affine_map, transform, change_of_basis, is_basis, translation_matrix,
-							  		 is_rotation_matrix, rotation_matrix, translation_rotation_matrix,
-							  		 change_basis_mesh, minimum_image_distance)
+from .quantities import electron_wavelength, interaction_parameter, lorentz
 from .voigt import gaussian, lorentzian, pseudo_voigt

skued/array_utils.py

@@ -3,7 +3,7 @@
 """
 
 from itertools import repeat
-from numpy import concatenate
+from numpy import concatenate, swapaxes
 
 def repeated_array(arr, num, axes = -1):
     """
@@ -44,4 +44,31 @@ def repeated_array(arr, num, axes = -1):
         for n, ax in zip(num[1:], axes[1:]):
             composite = concatenate(tuple(repeat(composite, times = n)), axis = ax)
 
-    return composite
+    return composite
+
+def mirror(arr, axes = None):
+    """ 
+    Reverse array over many axes. Generalization of arr[::-1] for many dimensions.
+
+    Parameters
+    ----------
+    arr : `~numpy.ndarray`
+        Array to be reversed
+    axes : int or tuple or None, optional
+        Axes to be reversed. Default is to reverse all axes.
+    
+    Returns
+    -------
+    out : 
+    """
+    if axes is None:
+        axes = range(arr.ndim)
+
+    elif isinstance(axes, int):
+        axes = (axes,)
+
+    # arr[::-1] reverses in the first axis direction
+    for axis in axes:
+        arr = swapaxes(swapaxes(arr, 0, axis)[::-1], 0, axis)
+
+    return arr

skued/quantities.py

@@ -4,10 +4,10 @@
 """
 import numpy as np
 
-c = 299792458           #m/s
-h = 6.63*10**(-34)      #J*s.
-e = 1.602*10**(-19)     #in C
-m0 = 9.109*10**(-31)    #in kg
+c = 299792458           # speed of light [m/s]
+h = 6.63*10**(-34)      # Planck's constant [J*s]
+e = 1.602*10**(-19)     # electron charge [C]
+m0 = 9.109*10**(-31)    # electron mass [kg]
 
 def lorentz(kV):
     """
@@ -16,7 +16,7 @@ def lorentz(kV):
     Parameters
     ----------
     kV : array_like or float
-        Electron gun voltage [kV].
+        Electron energy [kV].
     
     Returns
     -------
@@ -30,8 +30,8 @@ def electron_wavelength(kV):
         
     Parameters
     ----------
-    kV : float or ndarray
-        Voltage in kilovolts of the instrument
+    kV : array_like or float
+        Electron energy [kV].
     
     Returns
     -------
@@ -46,8 +46,8 @@ def interaction_parameter(kV):
 
     Parameters
     ----------
-    kV : float
-        Electron gun voltage [kV].
+    kV : array_like or float
+        Electron energy [kV].
     
     Returns
     -------
@@ -58,8 +58,6 @@ def interaction_parameter(kV):
     ----------
     .. Kirkland 2010 Eq. 5.6
     """
-    # Notes : check that interaction_parameter(100) == 0.92 * 1e-3
-    # See Kirkland 2010 p.84 for a discussion
     l = electron_wavelength(kV)
     V = kV * 1e3
 

skued/simulation/powdersim.py

@@ -5,7 +5,7 @@
 from ..voigt import pseudo_voigt
 import numpy as np
 
-def powdersim(crystal, scattering_length, broadening = True, **kwargs):
+def powdersim(crystal, scattering_length, fwhm_g = 0.01, fwhm_l = 0.02, **kwargs):
 	"""
 	Simulates polycrystalline diffraction pattern.
 
@@ -15,9 +15,10 @@ def powdersim(crystal, scattering_length, broadening = True, **kwargs):
 		Crystal from which to diffract.
 	scattering_length : `~numpy.ndarray`, shape (N,)
 		Range of scattering length over which to compute the diffraction pattern [2pi/Angs].
-	broadening: bool, optional
-		If True (default), returned pattern will display Voigt-related broadening.
-	
+	fwhm_g, fwhm_l : float, optional
+			Full-width at half-max of the Gaussian and Lorentzian parts of the Voigt profile.
+			See `skued.pseudo_voigt` for more details.
+
 	Returns
 	-------
 	pattern : `~numpy.ndarray`, shape (N,)
@@ -27,12 +28,10 @@ def powdersim(crystal, scattering_length, broadening = True, **kwargs):
 	scatt_length = np.sqrt(Gx**2 + Gy**2 + Gz**2)/(4*np.pi)
 	intensities = np.absolute(crystal.structure_factor((Gx, Gy, Gz)))**2
 
+	psf = pseudo_voigt(scattering_length, center = np.mean(scattering_length), fwhm_g = fwhm_g, fwhm_l = fwhm_l)
+
 	pattern = np.zeros_like(scattering_length)
-	if broadening:
-		for s, I in zip(scatt_length, intensities):
-			pattern += I * pseudo_voigt(scattering_length, s, 0.01, 0.02)
-	else:
-		for s, I in zip(scatt_length, intensities):
-			pattern += I * pseudo_voigt(scattering_length, s, 1e-5, 1e-5)
+	for s, i in zip(scatt_length, intensities):
+		pattern += i * pseudo_voigt(scattering_length, s, fwhm_g, fwhm_l)
 
 	return pattern

skued/structure/atom.py

@@ -91,8 +91,11 @@ def __init__(self, element, coords, displacement = None, **kwargs):
 
 		# Atomic potential parameters loaded on instantiation
 		# These are used to compute atomic potential
-		# TODO: add ref Kirkland 2010
-		_, a1, b1, a2, b2, a3, b3, c1, d1, c2, d2, c3, d3 = scattering_params[self.atomic_number]
+		try:
+			_, a1, b1, a2, b2, a3, b3, c1, d1, c2, d2, c3, d3 = scattering_params[self.atomic_number]
+		except KeyError:
+			raise ValueError('Scattering information for element {} is unavailable.'.format(self.element))
+
 		self._a = np.array((a1, a2, a3)).reshape((1,3))
 		self._b = np.array((b1, b2, b3)).reshape((1,3))
 		self._c = np.array((c1, c2, c3)).reshape((1,3))

skued/structure/cif_parser.py

@@ -85,18 +85,19 @@ class CIFParser(object):
 	"""
 	def __init__(self, filename, **kwargs):
 		""" 
+		Keyword arguments are passed to PyCIFRW's ReadCif function.
 		Parameters
 		----------
 		filename : str or path-like
 		"""
-		self.handle = open(filename, mode = 'r')
-		self.file = ReadCif(self.handle)
+		self._handle = open(filename, mode = 'r')
+		self.file = ReadCif(self._handle, **kwargs)
 
 	def __enter__(self):
 		return self
 
 	def __exit__(self, type, value, traceback):
-		self.handle.close()
+		self._handle.close()
 
 	@property
 	def _first_block(self):
@@ -107,9 +108,7 @@ def hall_symbol(self):
 		""" Returns the Hall symbol """
 		block = self._first_block
 
-		hall_symbol = None
-		for tag in ['_symmetry_space_group_name_Hall','_space_group_name_Hall']:
-			hall_symbol = block.get(tag) or hall_symbol
+		hall_symbol = block.get('_symmetry_space_group_name_Hall') or block.get('_space_group_name_Hall')
 
 		# In some rare cases, the given hall symbol in the file isn't standard,
 		# otherwise it would be a key in SymOpsHall
@@ -129,7 +128,7 @@ def hall_symbol(self):
 				hall_symbol = Number2Hall[table_number]
 
 		if hall_symbol is None:
-			raise ParseError('Hall number could not be inferred')
+			raise ParseError('Hall symbol could not be inferred')
 
 		if hall_symbol[0] == "-":
 			hall_symbol = "-" + hall_symbol[1].upper() + hall_symbol[2:].lower()

skued/structure/tests/test_atom.py

@@ -1,9 +1,11 @@
 # -*- coding: utf-8 -*-
-from random import choice
+from random import choice, seed
 import numpy as np
 from .. import Atom, atomic_number
 import unittest
 
+seed(23)
+
 ELEMENTS = list(atomic_number.keys())
 
 class TestAtom(unittest.TestCase):

skued/tests/test_array_utils.py

@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 import unittest
 import numpy as np
-from .. import repeated_array
+from .. import repeated_array, mirror
 
 np.random.seed(23)
 
@@ -33,5 +33,26 @@ def test_multiple_axes(self):
 			expected_new_shape = (self.arr.shape[0]*3, self.arr.shape[1]*2)
 			self.assertEqual(composite.shape, expected_new_shape)
 
+class TestMirror(unittest.TestCase):
+
+	def test_1D(self):
+		""" Test mirror() on a 1D array """
+		arr = np.zeros( (16,), dtype = np.float)
+		arr[15] = 1
+		self.assertTrue(np.allclose(arr[::-1], mirror(arr)))
+
+	def test_2D_all_axes(self):
+		""" Test mirror() on a 2D array for all axes """
+		arr = np.zeros( (16,16), dtype = np.float)
+		arr[15, 3] = 1
+		self.assertTrue(np.allclose(arr[::-1, ::-1], mirror(arr)))
+
+	def test_2D_one_axis(self):
+		""" Test mirror() on a 2D array for one axis """
+		arr = np.zeros( (16,16), dtype = np.float)
+		arr[15, 3] = 1
+		self.assertTrue(np.allclose(arr[:, ::-1], mirror(arr, axes = 1)))
+		self.assertTrue(np.allclose(arr[::-1, :], mirror(arr, axes = 0)))
+
 if __name__ == '__main__':
 	unittest.main()

skued/tests/test_quantities.py

@@ -18,10 +18,8 @@ class TestElectronWavelength(unittest.TestCase):
 
     def test_trivial(self):
         """ Test that the electron wavelength at zero energy is zero """
-        pass
-
-    def test_vectorized(self):
-        pass
+        self.assertAlmostEqual(electron_wavelength(10), 0.122128679)
+        self.assertAlmostEqual(electron_wavelength(200), 0.02509626)
 
 class TestInteractionParameter(unittest.TestCase):
 

skued/voigt.py

@@ -2,10 +2,9 @@
 """
 Voigt and pseudo-voigt curves, as well as related Gaussian and Lorentzian functions
 """
-from collections import Iterable
+from functools import lru_cache
 import numpy as np
 from numpy import pi
-from scipy.signal import fftconvolve
 
 def gaussian(coordinates, center, fwhm = None, std = None):
 	"""
@@ -109,6 +108,7 @@ def lorentzian(coordinates, center, fwhm):
 	factor = 1/(dim*pi)
 	return factor*core
 
+@lru_cache(maxsize = 16)
 def _pseudo_voigt_mixing_factor(width_l, width_g):
     """
     Returns the proportion of Lorentzian for the computation of a pseudo-Voigt profile.