Merge pull request #8 from LaurentRDC/develop

scikit-ued version 0.4.1

README.rst

@@ -20,9 +20,6 @@ To install the latest development version from `Github <https://github.com/Laure
 
     $ python -m pip install git+git://github.com/LaurentRDC/scikit-ued.git
 
-
-If your current Python installation doesn't have pip available, try `get-pip.py <bootstrap.pypa.io>`_
-
 After installing scikit-ued you can use it like any other Python module as :code:`skued`.
 
 API Reference

dev-requirements.txt

@@ -1,7 +1,6 @@
 numpy>=1.11
-nose==1.3.7
 coverage==4.3.1
 codecov==2.0.5
-Sphinx>=1.5.1
+Sphinx>=1.5.5
 cython>=0.25
 matplotlib

setup.py

@@ -1,12 +1,10 @@
 # -*- coding: utf-8 -*-
 import os
 import re
-import glob
-from itertools import chain
 from setuptools import setup, find_packages
-from Cython.Build import cythonize
+#from Cython.Build import cythonize
 
-base_package = 'skued'
+BASE_PACKAGE = 'skued'
 
 base_path = os.path.dirname(__file__)
 with open(os.path.join(base_path, 'skued', '__init__.py')) as f:
@@ -24,13 +22,9 @@
 with open('requirements.txt') as f:
     requirements = [line for line in f.read().split('\n') if len(line.strip())]
 
-packages = [base_package + '.' + x for x in find_packages(os.path.join(base_path, base_package))]
-if base_package not in packages:
-    packages.append(base_package)
-
-wavelets = chain.from_iterable([glob.glob('skued\\baseline\\data\\*.npy'), 
-                                glob.glob('skued\\baseline\\data\\*.npz')])
-
+packages = [BASE_PACKAGE + '.' + x for x in find_packages(os.path.join(base_path, BASE_PACKAGE))]
+if BASE_PACKAGE not in packages:
+    packages.append(BASE_PACKAGE)
 
 if __name__ == '__main__':
     setup(
@@ -39,6 +33,7 @@
         long_description = '\n\n'.join([readme, changes]),
         license = LICENSE,
         url = 'http://scikit-ued.readthedocs.io',
+        download_url = 'http://github.com/LaurentRDC/scikit-ued',
         version = VERSION,
         author = 'Laurent P. René de Cotret',
         author_email = 'laurent.renedecotret@mail.mcgill.ca',
@@ -47,12 +42,15 @@
         install_requires = requirements,
         keywords = ['skued'],
         packages = packages,
-		data_files = [('skued\\baseline\\data', wavelets)],
+        include_package_data = True,
         zip_safe = False,
-        classifiers = ['Intended Audience :: Developers',
+        classifiers = ['Environment :: Console',
+                       'Intended Audience :: Science/Research',
+                       'Topic :: Scientific'
                        'License :: OSI Approved :: MIT License',
                        'Natural Language :: English',
                        'Operating System :: OS Independent',
+                       'Programming Language :: Python',
                        'Programming Language :: Python :: 3.5',
                        'Programming Language :: Python :: 3.6']
     )

skued/__init__.py

@@ -2,8 +2,9 @@
 __author__ = 'Laurent P. René de Cotret'
 __email__ = 'laurent.renedecotret@mail.mcgill.ca'
 __license__ = 'MIT'
-__version__ = '0.4'
+__version__ = '0.4.1'
 
+from .array_utils import repeated_array
 from .parallel import pmap, preduce
 from .plot_utils import spectrum_colors
 from .quantities import lorentz, electron_wavelength, interaction_parameter

skued/array_utils.py

@@ -0,0 +1,47 @@
+"""
+Array utility functions
+"""
+
+from itertools import repeat
+from numpy import concatenate
+
+def repeated_array(arr, num, axes = -1):
+    """
+    Create a composite array from repeated copies of an array
+    
+    Parameters
+    ----------
+    arr : ndarray
+    
+    num : int or iterable of ints
+        Number of copies per axis. If provided as tuple, must be the same length
+        as 'axes' parameter. In case of `num` being 0 or an empty iterable, 
+        the inpur `arr` is returned.
+    axes : int or iterable of ints
+        Axis/axes over which to copy.
+    
+    Returns
+    -------
+    out : ndarray
+    
+    Raises
+    ------
+    ValueError
+        If num and axes are tuples of different lengths.
+    """
+    if not num:
+        return arr
+
+    if isinstance(num, int): num = (num,) 
+    if isinstance(axes, int): axes = (axes,)
+
+    if len(num) != len(axes):
+        raise ValueError('num and axes must have the same length')
+
+    composite = concatenate(tuple(repeat(arr, times = num[0])), axis = axes[0])
+
+    if len(num) > 1:
+        for n, ax in zip(num[1:], axes[1:]):
+            composite = concatenate(tuple(repeat(composite, times = n)), axis = ax)
+
+    return composite

skued/baseline/algorithms.py

@@ -58,7 +58,7 @@ def _iterative_baseline(array, max_iter, mask, background_regions, axes, approx_
 			background = np.swapaxes(np.swapaxes(background, 0, axis)[:-1], 0, axis)
 	return background
 
-def _dt_approx_rec(array, first_stage, wavelet, level, axis = -1):
+def _dt_approx_rec(array, first_stage, wavelet, mode, level, axis = -1):
     """
     Approximate reconstruction of a signal/image using the dual-tree approach.
     
@@ -74,6 +74,8 @@ def _dt_approx_rec(array, first_stage, wavelet, level, axis = -1):
         First-stage wavelet to use.
     wavelet : str, optional
         Complex wavelet to use in late stages.
+    mode : str, optional
+        Signal extension mode, see pywt.Modes.
     axis : int, optional
         Axis over which to compute the transform. Default is -1.
             
@@ -82,14 +84,14 @@ def _dt_approx_rec(array, first_stage, wavelet, level, axis = -1):
     reconstructed : ndarray
         Approximated reconstruction of the input array.
     """
-    coeffs = dtcwt(data = array, first_stage = first_stage, wavelet = wavelet, level = level, mode = 'constant', axis = axis)
+    coeffs = dtcwt(data = array, first_stage = first_stage, wavelet = wavelet, level = level, mode = mode, axis = axis)
     app_coeffs, *det_coeffs = coeffs
 
     det_coeffs = [np.zeros_like(det, dtype = np.complex) for det in det_coeffs]
     reconstructed = idtcwt(coeffs = [app_coeffs] + det_coeffs, first_stage = first_stage, wavelet = wavelet, mode = 'constant', axis = axis)
     return reconstructed
 
-def _dwt_approx_rec(array, level, wavelet, axis):
+def _dwt_approx_rec(array, level, wavelet, mode, axis):
 	"""
 	Approximate reconstruction of a signal/image. Uses the multi-level discrete wavelet 
 	transform to decompose a signal or an image, and reconstruct it using approximate 
@@ -100,13 +102,15 @@ def _dwt_approx_rec(array, level, wavelet, axis):
 	array : array-like
 		Array to be decomposed. Currently, only 1D and 2D arrays are supported.
 		Only even-lengths signals long the axis.
-	level : int or 'max' or None (deprecated)
+	level : int or None (deprecated)
 		Decomposition level. A higher level will result in a coarser approximation of
 		the input array. If the level is higher than the maximum possible decomposition level,
 		the maximum level is used.
 		If None, the maximum possible decomposition level is used.
 	wavelet : str or Wavelet object
 		Can be any argument accepted by PyWavelet.Wavelet, e.g. 'db10'
+    mode : str, optional
+        Signal extension mode, see pywt.Modes.
     axis : int, optional
         Axis over which to compute the transform. Default is -1.
             
@@ -135,15 +139,15 @@ def _dwt_approx_rec(array, level, wavelet, axis):
 	# For 2D array, check the condition with shortest dimension min(array.shape). This is how
 	# it is done in PyWavelet.wavedec2.
 	max_level = pywt.dwt_max_level(data_len = array.shape[axis], filter_len = wavelet.dec_len)
-	if level is None or level is 'max':
+	if level is None:
 		level = max_level
 	elif max_level < level:
 		warn('Decomposition level {} higher than maximum {}. Maximum is used.'.format(level, max_level))
 		level = max_level
 
 	# By now, we are sure that the decomposition level will be supported.
 	# Decompose the signal using the multilevel discrete wavelet transform
-	coeffs = pywt.wavedec(data = array, wavelet = wavelet, level = level, mode = 'constant', axis = axis)
+	coeffs = pywt.wavedec(data = array, wavelet = wavelet, level = level, mode = mode, axis = axis)
 	app_coeffs, det_coeffs = coeffs[0], coeffs[1:]
 
 	# Replace detail coefficients by 0; keep the correct length so that the
@@ -157,7 +161,7 @@ def _dwt_approx_rec(array, level, wavelet, axis):
 		reconstructed = np.swapaxes(np.swapaxes(reconstructed, 0, axis)[:array.shape[axis]], 0, axis)
 	return  reconstructed
 
-def _dwt_approx_rec2(array, level, wavelet, axis):
+def _dwt_approx_rec2(array, level, wavelet, mode, axis):
 	"""
 	Approximate reconstruction of a signal/image. Uses the multi-level discrete wavelet 
 	transform to decompose a signal or an image, and reconstruct it using approximate 
@@ -175,6 +179,8 @@ def _dwt_approx_rec2(array, level, wavelet, axis):
 		If None, the maximum possible decomposition level is used.
 	wavelet : str or Wavelet object
 		Can be any argument accepted by PyWavelet.Wavelet, e.g. 'db10'
+    mode : str, optional
+        Signal extension mode, see pywt.Modes.
 	axis : 2-tuple of ints
             
 	Returns
@@ -208,15 +214,15 @@ def _dwt_approx_rec2(array, level, wavelet, axis):
 
 	# By now, we are sure that the decomposition level will be supported.
 	# Decompose the signal using the multilevel discrete wavelet transform
-	coeffs = pywt.wavedec2(data = array, wavelet = wavelet, level = level, mode = 'constant', axes = axis)
+	coeffs = pywt.wavedec2(data = array, wavelet = wavelet, level = level, mode = mode, axes = axis)
 	app_coeffs, det_coeffs = coeffs[0], coeffs[1:]
 
 	# Replace detail coefficients by 0; keep the correct length so that the
 	# reconstructed signal has the same size as the (possibly upsampled) signal
 	# The structure of coefficients depends on the dimensionality	
 	# Reconstruct signal
 	reconstructed = pywt.waverec2([app_coeffs] + [(None, None, None)]*len(det_coeffs), 
-								  wavelet = wavelet, mode = 'constant', axes = axis)
+								  wavelet = wavelet, mode = mode, axes = axis)
 
 	# Sometimes pywt.waverec returns a signal that is longer than the original signal
 	for ax in axis:
@@ -225,7 +231,7 @@ def _dwt_approx_rec2(array, level, wavelet, axis):
 	return reconstructed
 
 def baseline_dt(array, max_iter, level = None, first_stage = 'sym6', wavelet = 'qshift1', 
-				background_regions = [], mask = None, axis = -1):
+				background_regions = [], mask = None, mode = 'constant', axis = -1):
 	"""
 	Iterative method of baseline determination based on the dual-tree complex wavelet transform.
 	This function only works in 1D, along an axis. For baseline of 2D arrays, see baseline_dwt.
@@ -236,9 +242,9 @@ def baseline_dt(array, max_iter, level = None, first_stage = 'sym6', wavelet = '
 		Data with background. Can be either 1D signal or 2D array.
 	max_iter : int
 		Number of iterations to perform.
-	level : int or 'max', optional
+	level : int or None, optional
 		Decomposition level. A higher level will result in a coarser approximation of
-		the input signal (read: a lower frequency baseline). If 'max' (default), the maximum level
+		the input signal (read: a lower frequency baseline). If None (default), the maximum level
 		possible is used.
 	first_stage : str, optional
 		Wavelet to use for the first stage. See skued.baseline.ALL_FIRST_STAGE for a list of suitable arguments
@@ -261,6 +267,8 @@ def baseline_dt(array, max_iter, level = None, first_stage = 'sym6', wavelet = '
     
 	mask : ndarray, dtype bool, optional
 		Mask array that evaluates to True for pixels that are invalid. 
+    mode : str, optional
+        Signal extension mode, see pywt.Modes.
 	axis : int, optional
 		Axis over which to compute the wavelet transform. Default is -1
     
@@ -276,10 +284,11 @@ def baseline_dt(array, max_iter, level = None, first_stage = 'sym6', wavelet = '
 	"""
 	return _iterative_baseline(array = array, max_iter = max_iter, background_regions = background_regions,
 									mask = mask, axes = (axis,), approx_rec_func = _dt_approx_rec,
-									func_kwargs = {'level': level, 'wavelet': wavelet,
+									func_kwargs = {'level': level, 'wavelet': wavelet, 'mode': mode,
 													'first_stage': first_stage, 'axis': axis})
 
-def baseline_dwt(array, max_iter, level = None, wavelet = 'sym6', background_regions = [], mask = None, axis = -1):
+def baseline_dwt(array, max_iter, level = None, wavelet = 'sym6', background_regions = [], 
+				 mask = None, mode = 'constant', axis = -1):
 	"""
 	Iterative method of baseline determination, based on the discrete wavelet transform. 
     
@@ -313,6 +322,8 @@ def baseline_dwt(array, max_iter, level = None, wavelet = 'sym6', background_reg
 	mask : ndarray, dtype bool, optional
 		Mask array that evaluates to True for pixels that are invalid. Useful to determine which pixels are masked
 		by a beam block.
+    mode : str, optional
+        Signal extension mode, see pywt.Modes.
 	axis : int or tuple, optional
 		Axis over which to compute the wavelet transform. Can also be a 2-tuple of ints for 2D baseline
     
@@ -333,4 +344,4 @@ def baseline_dwt(array, max_iter, level = None, wavelet = 'sym6', background_reg
 
 	return _iterative_baseline(array, max_iter = max_iter, background_regions = background_regions, 
 							   mask = mask, axes = axis, approx_rec_func = approx_rec_func[len(axis)], 
-							   func_kwargs = {'level': level, 'wavelet': wavelet, 'axis': axis})
+							   func_kwargs = {'level': level, 'wavelet': wavelet, 'axis': axis, 'mode':mode})

skued/baseline/tests/test_algorithms.py

@@ -19,28 +19,28 @@ def test_approx_rec(self):
 
 		with self.subTest('1D shape'):
 			arr = np.random.random(size = (102,))
-			rec_arr = _dwt_approx_rec(arr, level = 2, wavelet = 'db6', axis = -1)
+			rec_arr = _dwt_approx_rec(arr, level = 2, wavelet = 'db6', mode = 'constant', axis = -1)
 			self.assertSequenceEqual(rec_arr.shape, arr.shape)
 
 		with self.subTest('1D along axis'):
 			arr2 = np.random.random(size = (21, 52))
-			rec_arr2 = _dwt_approx_rec(arr2, level = 2, wavelet = 'sym4', axis = 1)
+			rec_arr2 = _dwt_approx_rec(arr2, level = 2, wavelet = 'sym4', mode = 'constant', axis = 1)
 			self.assertSequenceEqual(rec_arr2.shape, arr2.shape)
 
 		with self.subTest('2D shape'):
 			arr2 = np.random.random(size = (22,52))
-			rec_arr2 = _dwt_approx_rec2(arr2, level = 2, wavelet = 'sym6', axis = (-2, -1))
+			rec_arr2 = _dwt_approx_rec2(arr2, level = 2, wavelet = 'sym6', mode = 'constant', axis = (-2, -1))
 			self.assertSequenceEqual(rec_arr2.shape, arr2.shape)
 
 		with self.subTest('2D along axes'):
 			arr2 = np.random.random(size = (22,52, 5))
-			rec_arr2 = _dwt_approx_rec2(arr2, level = 2, wavelet = 'sym6', axis = (0, 1))
+			rec_arr2 = _dwt_approx_rec2(arr2, level = 2, wavelet = 'sym6', mode = 'constant', axis = (0, 1))
 			self.assertSequenceEqual(rec_arr2.shape, arr2.shape)
 
 	def test_trivial_case_1d(self):
 		""" The baseline for a 1d array of zeros should be zeros """
 		arr = np.zeros_like(self.arr)
-		self.assertTrue(np.allclose(arr, baseline_dwt(arr, max_iter = 10, level = 'max')))
+		self.assertTrue(np.allclose(arr, baseline_dwt(arr, max_iter = 10, level = None)))
 
 	def test_trivial_case_2d(self):
 		""" The baseline for a 2d array of zeros should be zeros """
@@ -76,12 +76,12 @@ def test_approx_rec(self):
 
 		with self.subTest('1D shape'):
 			arr = np.random.random(size = (102,))
-			rec_arr = _dt_approx_rec(arr, level = 2, first_stage = 'db1', wavelet = 'qshift3', axis = -1)
+			rec_arr = _dt_approx_rec(arr, level = 2, first_stage = 'db1', wavelet = 'qshift3', mode = 'smooth', axis = -1)
 			self.assertSequenceEqual(rec_arr.shape, arr.shape)
 
 		with self.subTest('2D along axis'):
 			arr2 = np.random.random(size = (21, 52))
-			rec_arr2 = _dt_approx_rec(arr2, level = 2, first_stage = 'db1', wavelet = 'qshift3', axis = 1)
+			rec_arr2 = _dt_approx_rec(arr2, level = 2, first_stage = 'db1', wavelet = 'qshift3', mode = 'smooth', axis = 1)
 			self.assertSequenceEqual(rec_arr2.shape, arr2.shape)
 
 	def test_trivial_case(self):