# pbrod/scipy forked from scipy/scipy

Merge branch 'master' of https://github.com/scipy/scipy into ticket1842

2 parents 95ca499 + d76b18c commit 01b7092fd9869cc95e0e1c8a63d8330fdd3fb1e3 committed Mar 16, 2013
 @@ -29,6 +29,27 @@ Deprecated features Backwards incompatible changes ============================== +LIL matrix assignment +--------------------- +Assigning values to LIL matrices with two index arrays now works similarly as +assigning into ndarrays:: + + >>> x = lil_matrix((3, 3)) + >>> x[[0,1,2],[0,1,2]]=[0,1,2] + >>> x.todense() + matrix([[ 0., 0., 0.], + [ 0., 1., 0.], + [ 0., 0., 2.]]) + +rather than giving the result:: + + >>> x.todense() + matrix([[ 0., 1., 2.], + [ 0., 1., 2.], + [ 0., 1., 2.]]) + +Users relying on the previous behavior will need to revisit their code. + Other changes =============
 @@ -109,12 +109,16 @@ class UnivariateSpline(object): -------- >>> from numpy import linspace,exp >>> from numpy.random import randn + >>> import matplotlib.pyplot as plt >>> from scipy.interpolate import UnivariateSpline >>> x = linspace(-3, 3, 100) >>> y = exp(-x**2) + randn(100)/10 >>> s = UnivariateSpline(x, y, s=1) >>> xs = linspace(-3, 3, 1000) >>> ys = s(xs) + >>> plt.plot(x, y, '.-') + >>> plt.plot(xs, ys) + >>> plt.show() xs,ys is now a smoothed, super-sampled version of the noisy gaussian x,y. @@ -315,11 +319,15 @@ class InterpolatedUnivariateSpline(UnivariateSpline): >>> from numpy import linspace,exp >>> from numpy.random import randn >>> from scipy.interpolate import InterpolatedUnivariateSpline + >>> import matplotlib.pyplot as plt >>> x = linspace(-3, 3, 100) >>> y = exp(-x**2) + randn(100)/10 >>> s = InterpolatedUnivariateSpline(x, y) >>> xs = linspace(-3, 3, 1000) >>> ys = s(xs) + >>> plt.plot(x, y, '.-') + >>> plt.plot(xs, ys) + >>> plt.show() xs,ys is now a smoothed, super-sampled version of the noisy gaussian x,y @@ -357,7 +365,7 @@ class LSQUnivariateSpline(UnivariateSpline): Input dimension of data points -- must be increasing y : (N,) array_like Input dimension of data points - t: (M,) array_like + t : (M,) array_like interior knots of the spline. Must be in ascending order and bbox[0]>> from numpy import linspace,exp >>> from numpy.random import randn >>> from scipy.interpolate import LSQUnivariateSpline + >>> import matplotlib.pyplot as plt >>> x = linspace(-3,3,100) >>> y = exp(-x**2) + randn(100)/10 >>> t = [-1,0,1] >>> s = LSQUnivariateSpline(x,y,t) >>> xs = linspace(-3,3,1000) >>> ys = s(xs) - + >>> plt.plot(x, y, '.-') + >>> plt.plot(xs, ys) + >>> plt.show() + xs,ys is now a smoothed, super-sampled version of the noisy gaussian x,y with knots [-3,-1,0,1,3]
 @@ -520,6 +520,14 @@ def write_element(self, arr, mdtype=None): ''' write tag and data ''' if mdtype is None: mdtype = NP_TO_MTYPES[arr.dtype.str[1:]] + + # We are writing a little-endian Matlab file but our incoming arrays may + # be big-endian. In particular, they might be big-endian because we originally + # *read* them from a big-endian Matlab file + byte_order = arr.dtype.byteorder + if byte_order == '>' or (byte_order == '=' and not boc.sys_is_le): + arr = arr.byteswap().newbyteorder() + byte_count = arr.size*arr.itemsize if byte_count <= 4: self.write_smalldata_element(arr, mdtype, byte_count)
 @@ -203,11 +203,10 @@ cdef class VarReader5: if isinstance(key, str): continue self.class_dtypes[key] = dt - # cache correctly byte ordered dtypes - if self.little_endian: - self.U1_dtype = np.dtype('U1') + # Always use U1 rather than U1 for interpreting string + # data because the strings are created by the Python runtime + # by .decode() and hence use native byte order rather the mat file's + self.U1_dtype = np.dtype('U1') bool_dtype = np.dtype('bool') def set_stream(self, fobj):
Binary file not shown.
 @@ -818,6 +818,34 @@ def test_empty_string(): stream.close() +def test_read_big_endian(): + # make sure big-endian data is read correctly + estring_fname = pjoin(test_data_path, 'big_endian.mat') + fp = open(estring_fname, 'rb') + rdr = MatFile5Reader_future(fp) + d = rdr.get_variables() + fp.close() + assert_array_equal(d['strings'], np.array([[u'hello'], + [u'world']], dtype=np.object)) + assert_array_equal(d['floats'], np.array([[ 2., 3.], + [ 3., 4.]], dtype=np.float32)) + + +def test_write_big_endian(): + # we don't support writing actual big-endian .mat files, but we need to + # behave correctly if the user supplies a big-endian numpy array to write out + stream = BytesIO() + savemat_future(stream, {'a': np.array([[ 2., 3.], + [ 3., 4.]], dtype='>f4'), + 'b': np.array([u'hello', u'world'], dtype='>U')}) + rdr = MatFile5Reader_future(stream) + d = rdr.get_variables() + assert_array_equal(d['a'], np.array([[ 2., 3.], + [ 3., 4.]], dtype='f4')) + assert_array_equal(d['b'], np.array([u'hello', u'world'], dtype='U')) + stream.close() + + def test_mat4_3d(): # test behavior when writing 3D arrays to matlab 4 files stream = BytesIO()
 @@ -4,7 +4,7 @@ from numpy import array -def imread(fname, flatten=False): +def imread(fname, flatten=False, mode=None): """ Load an image from file. @@ -14,6 +14,9 @@ def imread(fname, flatten=False): Image file name, e.g. ``test.jpg``. flatten : bool, optional If true, convert the output to grey-scale. Default is False. + mode : str, optional + mode to convert image to, e.g. ``RGB``. + Returns ------- @@ -36,10 +39,11 @@ def imread(fname, flatten=False): " http://pypi.python.org/pypi/PIL/ for installation" " instructions.") - fp = open(fname, "rb") - im = Image.open(fp) - if flatten: - im = im.convert('F') - result = array(im) - fp.close() - return result + with open(fname, "rb") as fp: + im = Image.open(fp) + if mode: + im = im.convert(mode) + if flatten: + im = im.convert('F') + result = array(im) + return result
 @@ -14,7 +14,7 @@ @dec.skipif(pil_missing, msg="The Python Image Library could not be found.") def test_imread(): lp = os.path.join(os.path.dirname(__file__), 'dots.png') - img = ndi.imread(lp) + img = ndi.imread(lp, mode="RGB") assert_array_equal(img.shape, (300, 420, 3)) img = ndi.imread(lp, flatten=True)
 @@ -28,6 +28,34 @@ class BadCoefficients(UserWarning): def findfreqs(num, den, N): + """ + Find an array of frequencies for computing the response of a filter. + + Parameters + ---------- + num, den : array_like, 1-D + The polynomial coefficients of the numerator and denominator of the + transfer function of the filter or LTI system. The coefficients are + ordered from highest to lowest degree. + N : int + The length of the array to be computed. + + Returns + ------- + w : (N,) ndarray + A 1-D array of frequencies, logarithmically spaced. + + Examples + -------- + Find a set of nine frequencies that span the "interesting part" of the + frequency response for the filter with the transfer function + H(s) = s / (s^2 + 8s + 25) + + >>> findfreqs([1, 0], [1, 8, 25], N=9) + array([ 1.00000000e-02, 3.16227766e-02, 1.00000000e-01, + 3.16227766e-01, 1.00000000e+00, 3.16227766e+00, + 1.00000000e+01, 3.16227766e+01, 1.00000000e+02]) + """ ep = atleast_1d(roots(den)) + 0j tz = atleast_1d(roots(num)) + 0j
 @@ -613,44 +613,6 @@ def _default_response_times(A, n): return t -def _default_response_frequencies(A, n): - """ - Compute a reasonable set of frequency points for Bode plot. - - This function is used by `bode` to compute the frequency points (in rad/s) - when the `w` argument to the function is None. - - Parameters - ---------- - A : ndarray - The system matrix, which is square. - n : int - The number of frequency samples to generate. - - Returns - ------- - w : ndarray - The 1-D array of length `n` of frequency samples (in rad/s) at which - the response is to be computed. - """ - vals = linalg.eigvals(A) - # Remove poles at 0 because they don't help us determine an interesting - # frequency range. (And if we pass a 0 to log10() below we will crash.) - poles = [pole for pole in vals if pole != 0] - # If there are no non-zero poles, just hardcode something. - if len(poles) == 0: - minpole = 1 - maxpole = 1 - else: - minpole = min(abs(real(poles))) - maxpole = max(abs(real(poles))) - # A reasonable frequency range is two orders of magnitude before the - # minimum pole (slowest) and two orders of magnitude after the maximum pole - # (fastest). - w = numpy.logspace(numpy.log10(minpole) - 2, numpy.log10(maxpole) + 2, n) - return w - - def impulse(system, X0=None, T=None, N=None): """Impulse response of continuous-time system. @@ -936,14 +898,14 @@ def bode(system, w=None, n=100): sys = lti(*system) if w is None: - w = _default_response_frequencies(sys.A, n) + worN = n else: - w = numpy.asarray(w) + worN = w + w, y = freqs(sys.num, sys.den, worN=worN) - jw = w * 1j - y = numpy.polyval(sys.num, jw) / numpy.polyval(sys.den, jw) mag = 20.0 * numpy.log10(abs(y)) phase = numpy.arctan2(y.imag, y.real) * 180.0 / numpy.pi + return w, mag, phase
 @@ -725,6 +725,30 @@ def lfiltic(b, a, y, x=None): def deconvolve(signal, divisor): """Deconvolves `divisor` out of `signal`. + Parameters + ---------- + signal : array + Signal input + divisor : array + Divisor input + + Returns + ------- + q : array + Quotient of the division + r : array + Remainder + + Examples + -------- + >>> from scipy import signal + >>> sig = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1,]) + >>> filter = np.array([1,1,0]) + >>> res = signal.convolve(sig, filter) + >>> signal.deconvolve(res, filter) + (array([ 0., 0., 0., 0., 0., 1., 1., 1., 1.]), + array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])) + """ num = atleast_1d(signal) den = atleast_1d(divisor)
 @@ -285,19 +285,15 @@ def test_04(self): assert_almost_equal(phase, expected_phase) def test_05(self): - """Test that bode() finds a reasonable frequency range. - - A reasonable frequency range is two orders of magnitude before the - minimum (slowest) pole and two orders of magnitude after the maximum - (fastest) pole. - """ + """Test that bode() finds a reasonable frequency range.""" # 1st order low-pass filter: H(s) = 1 / (s + 1) system = lti([1], [1, 1]) vals = linalg.eigvals(system.A) minpole = min(abs(np.real(vals))) maxpole = max(abs(np.real(vals))) n = 10; - expected_w = np.logspace(np.log10(minpole) - 2, np.log10(maxpole) + 2, n) + # Expected range is from 0.01 to 10. + expected_w = np.logspace(-2, 1, n) w, mag, phase = bode(system, n=n) assert_almost_equal(w, expected_w) @@ -308,6 +304,12 @@ def test_06(self): w, mag, phase = bode(system, n=2) assert_equal(w[0], 0.01) # a fail would give not-a-number + def test_07(self): + """bode() should not fail on a system with pure imaginary poles.""" + # The test passes if bode doesn't raise an exception. + system = lti([1], [1, 0, 100]) + w, mag, phase = bode(system, n=2) + class Test_freqresp(object): @@ -377,5 +379,6 @@ def test_pole_zero(self): w, H = freqresp(system, n=2) assert_equal(w[0], 0.01) # a fail would give not-a-number + if __name__ == "__main__": run_module_suite()
 @@ -10,15 +10,15 @@ import numpy as np -from .data import _data_matrix +from .data import _data_matrix, _minmax_mixin from .compressed import _cs_matrix from .base import isspmatrix, _formats from .sputils import isshape, getdtype, to_native, upcast from . import sparsetools from .sparsetools import bsr_matvec, bsr_matvecs, csr_matmat_pass1, \ bsr_matmat_pass2, bsr_transpose, bsr_sort_indices -class bsr_matrix(_cs_matrix): +class bsr_matrix(_cs_matrix, _minmax_mixin): """Block Sparse Row matrix This can be instantiated in several ways: