def test_fix_with_subclass(self):
        class MyArray(nx.ndarray):
            def __new__(cls, data, metadata=None):
                res = nx.array(data, copy=True).view(cls)
                res.metadata = metadata
                return res

            def __array_wrap__(self, obj, context=None):
                if isinstance(obj, MyArray):
                    obj.metadata = self.metadata
                return obj

            def __array_finalize__(self, obj):
                self.metadata = getattr(obj, 'metadata', None)
                return self

        a = nx.array([1.1, -1.1])
        m = MyArray(a, metadata='foo')
        f = ufl.fix(m)
        assert_array_equal(f, nx.array([1, -1]))
        assert_(isinstance(f, MyArray))
        assert_equal(f.metadata, 'foo')

        # check 0d arrays don't decay to scalars
        m0d = m[0,...]
        m0d.metadata = 'bar'
        f0d = ufl.fix(m0d)
        assert_(isinstance(f0d, MyArray))
        assert_equal(f0d.metadata, 'bar')

def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                         header=''):
    from numpy.core import array, isnan, any
    x = array(x, copy=False, subok=True)
    y = array(y, copy=False, subok=True)

    def isnumber(x):
        return x.dtype.char in '?bhilqpBHILQPfdgFDG'

    try:
        cond = (x.shape==() or y.shape==()) or x.shape == y.shape
        if not cond:
            msg = build_err_msg([x, y],
                                err_msg
                                + '\n(shapes %s, %s mismatch)' % (x.shape,
                                                                  y.shape),
                                verbose=verbose, header=header,
                                names=('x', 'y'))
            if not cond :
                raise AssertionError(msg)

        if (isnumber(x) and isnumber(y)) and (any(isnan(x)) or any(isnan(y))):
            # Handling nan: we first check that x and y have the nan at the
            # same locations, and then we mask the nan and do the comparison as
            # usual.
            xnanid = isnan(x)
            ynanid = isnan(y)
            try:
                assert_array_equal(xnanid, ynanid)
            except AssertionError:
                msg = build_err_msg([x, y],
                                    err_msg
                                    + '\n(x and y nan location mismatch %s, ' \
                                    '%s mismatch)' % (xnanid, ynanid),
                                    verbose=verbose, header=header,
                                    names=('x', 'y'))
            val = comparison(x[~xnanid], y[~ynanid])
        else:
            val = comparison(x,y)
        if isinstance(val, bool):
            cond = val
            reduced = [0]
        else:
            reduced = val.ravel()
            cond = reduced.all()
            reduced = reduced.tolist()
        if not cond:
            match = 100-100.0*reduced.count(1)/len(reduced)
            msg = build_err_msg([x, y],
                                err_msg
                                + '\n(mismatch %s%%)' % (match,),
                                verbose=verbose, header=header,
                                names=('x', 'y'))
            if not cond :
                raise AssertionError(msg)
    except ValueError:
        msg = build_err_msg([x, y], err_msg, verbose=verbose, header=header,
                            names=('x', 'y'))
        raise ValueError(msg)

 def test_3D_array(self):
     a = array([[1, 2], [1, 2]])
     b = array([[2, 3], [2, 3]])
     a = array([a, a])
     b = array([b, b])
     res = map(atleast_1d, [a, b])
     desired = [a, b]
     assert_array_equal(res, desired)

 def test_3D_array(self):
     a = array([[1, 2], [1, 2]])
     b = array([[2, 3], [2, 3]])
     a = array([a, a])
     b = array([b, b])
     res = [atleast_2d(a), atleast_2d(b)]
     desired = [a, b]
     assert_array_equal(res, desired)

    def test_bad_out_shape(self):
        a = array([1, 2])
        b = array([3, 4])

        assert_raises(ValueError, concatenate, (a, b), out=np.empty(5))
        assert_raises(ValueError, concatenate, (a, b), out=np.empty((4,1)))
        assert_raises(ValueError, concatenate, (a, b), out=np.empty((1,4)))
        concatenate((a, b), out=np.empty(4))

 def test_log2(self):
     a = nx.array([4.5, 2.3, 6.5])
     out = nx.zeros(a.shape, float)
     tgt = nx.array([2.169925, 1.20163386, 2.70043972])
     res = ufl.log2(a)
     assert_almost_equal(res, tgt)
     res = ufl.log2(a, out)
     assert_almost_equal(res, tgt)
     assert_almost_equal(out, tgt)

def test_stack():
    # non-iterable input
    assert_raises(TypeError, stack, 1)

    # 0d input
    for input_ in [(1, 2, 3),
                   [np.int32(1), np.int32(2), np.int32(3)],
                   [np.array(1), np.array(2), np.array(3)]]:
        assert_array_equal(stack(input_), [1, 2, 3])
    # 1d input examples
    a = np.array([1, 2, 3])
    b = np.array([4, 5, 6])
    r1 = array([[1, 2, 3], [4, 5, 6]])
    assert_array_equal(np.stack((a, b)), r1)
    assert_array_equal(np.stack((a, b), axis=1), r1.T)
    # all input types
    assert_array_equal(np.stack(list([a, b])), r1)
    assert_array_equal(np.stack(array([a, b])), r1)
    # all shapes for 1d input
    arrays = [np.random.randn(3) for _ in range(10)]
    axes = [0, 1, -1, -2]
    expected_shapes = [(10, 3), (3, 10), (3, 10), (10, 3)]
    for axis, expected_shape in zip(axes, expected_shapes):
        assert_equal(np.stack(arrays, axis).shape, expected_shape)
    assert_raises_regex(np.AxisError, 'out of bounds', stack, arrays, axis=2)
    assert_raises_regex(np.AxisError, 'out of bounds', stack, arrays, axis=-3)
    # all shapes for 2d input
    arrays = [np.random.randn(3, 4) for _ in range(10)]
    axes = [0, 1, 2, -1, -2, -3]
    expected_shapes = [(10, 3, 4), (3, 10, 4), (3, 4, 10),
                       (3, 4, 10), (3, 10, 4), (10, 3, 4)]
    for axis, expected_shape in zip(axes, expected_shapes):
        assert_equal(np.stack(arrays, axis).shape, expected_shape)
    # empty arrays
    assert_(stack([[], [], []]).shape == (3, 0))
    assert_(stack([[], [], []], axis=1).shape == (0, 3))
    # out
    out = np.zeros_like(r1)
    np.stack((a, b), out=out)
    assert_array_equal(out, r1)
    # edge cases
    assert_raises_regex(ValueError, 'need at least one array', stack, [])
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [1, np.arange(3)])
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.arange(3), 1])
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.arange(3), 1], axis=1)
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.zeros((3, 3)), np.zeros(3)], axis=1)
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.arange(2), np.arange(3)])
    # generator is deprecated
    with assert_warns(FutureWarning):
        result = stack((x for x in range(3)))
    assert_array_equal(result, np.array([0, 1, 2]))

    def test_isneginf(self):
        a = nx.array([nx.inf, -nx.inf, nx.nan, 0.0, 3.0, -3.0])
        out = nx.zeros(a.shape, bool)
        tgt = nx.array([False, True, False, False, False, False])

        res = ufl.isneginf(a)
        assert_equal(res, tgt)
        res = ufl.isneginf(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)

    def test_fix(self):
        a = nx.array([[1.0, 1.1, 1.5, 1.8], [-1.0, -1.1, -1.5, -1.8]])
        out = nx.zeros(a.shape, float)
        tgt = nx.array([[1., 1., 1., 1.], [-1., -1., -1., -1.]])

        res = ufl.fix(a)
        assert_equal(res, tgt)
        res = ufl.fix(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)
        assert_equal(ufl.fix(3.14), 3)

def Trebuchet1(t, state, param):
	g = param
	Theta + dTheta = state
	f   = array([(A*dTheta**2*h*n - A*dTheta**2*h*s - A*g*n + A*g*s - 2*M*g*n + 2*g*m*s)*cos(Theta)/2])
	m   = array([
		[A*(n + s)**2/12 + A*(4*h**2 - 4*h*(n - s)*sin(Theta) + (n - s)**2)/4 + M*n**2 + m*s**2],
		])
	n   = array([0])
	if h - s*sin(Theta) < 0:
		n = array([-g*m*s*cos(Theta)])
	f = f + n
	return concatenate((solve(m,f),[dTheta]))

def test_stack():
    # non-iterable input
    assert_raises(TypeError, stack, 1)

    # 0d input
    for input_ in [(1, 2, 3),
                   [np.int32(1), np.int32(2), np.int32(3)],
                   [np.array(1), np.array(2), np.array(3)]]:
        assert_array_equal(stack(input_), [1, 2, 3])
    # 1d input examples
    a = np.array([1, 2, 3])
    b = np.array([4, 5, 6])
    r1 = array([[1, 2, 3], [4, 5, 6]])
    assert_array_equal(np.stack((a, b)), r1)
    assert_array_equal(np.stack((a, b), axis=1), r1.T)
    # all input types
    assert_array_equal(np.stack(list([a, b])), r1)
    assert_array_equal(np.stack(array([a, b])), r1)
    # all shapes for 1d input
    arrays = [np.random.randn(3) for _ in range(10)]
    axes = [0, 1, -1, -2]
    expected_shapes = [(10, 3), (3, 10), (3, 10), (10, 3)]
    for axis, expected_shape in zip(axes, expected_shapes):
        assert_equal(np.stack(arrays, axis).shape, expected_shape)
    assert_raises_regex(IndexError, 'out of bounds', stack, arrays, axis=2)
    assert_raises_regex(IndexError, 'out of bounds', stack, arrays, axis=-3)
    # all shapes for 2d input
    arrays = [np.random.randn(3, 4) for _ in range(10)]
    axes = [0, 1, 2, -1, -2, -3]
    expected_shapes = [(10, 3, 4), (3, 10, 4), (3, 4, 10),
                        (3, 4, 10), (3, 10, 4), (10, 3, 4)]
    for axis, expected_shape in zip(axes, expected_shapes):
        assert_equal(np.stack(arrays, axis).shape, expected_shape)
    # empty arrays
    assert_(stack([[], [], []]).shape == (3, 0))
    assert_(stack([[], [], []], axis=1).shape == (0, 3))
    # edge cases
    assert_raises_regex(ValueError, 'need at least one array', stack, [])
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [1, np.arange(3)])
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.arange(3), 1])
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.arange(3), 1], axis=1)
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.zeros((3, 3)), np.zeros(3)], axis=1)
    assert_raises_regex(ValueError, 'must have the same shape',
                        stack, [np.arange(2), np.arange(3)])
    # np.matrix
    m = np.matrix([[1, 2], [3, 4]])
    assert_raises_regex(ValueError, 'shape too large to be a matrix',
                        stack, [m, m])

def Trebuchet2(t, state, param):
	g = param
	Theta + dTheta = state
	f   = array([A*dTheta**2*h*(n - s)*cos(Theta)/2 - A*g*(n - s)*cos(Theta)/2 - M*dPhi*dTheta*n*p*cos(Phi) - M*dPhi*n*p*(dPhi + dTheta)*cos(Phi) - M*g*(n*cos(Theta) + p*(sin(Phi)*cos(Theta) + sin(Theta)*cos(Phi))) + g*m*s*cos(Theta), M*p*(dTheta**2*n*cos(Phi) - g*sin(Phi)*cos(Theta) - g*sin(Theta)*cos(Phi))])
	m   = array([
		[A*(n + s)**2/12 + A*(4*h**2 - 4*h*(n - s)*sin(Theta) + (n - s)**2)/4 + M*n**2 + 2*M*n*p*sin(Phi) + M*p**2 + m*s**2, M*p*(n*sin(Phi) + p)],
		[M*p*(n*sin(Phi) + p), M*p**2],
		])
	n   = array([0, 0])
	if h - s*sin(Theta) < 0:
		n = array([-g*m*s*cos(Theta), 0])
	f = f + n
	return concatenate((solve(m,f),[dTheta, dPhi]))

def Trebuchet3(t, state, param):
	g = param
	Theta + dTheta = state
	f   = array([A*dTheta**2*h*(n - s)*cos(Theta)/2 - A*g*(n - s)*cos(Theta)/2 - M*g*n*cos(Theta) + dPsi*dTheta*m*q*s*cos(Psi) + dPsi*m*q*s*(dPsi + dTheta)*cos(Psi) - g*m*(q*(sin(Psi)*cos(Theta) + sin(Theta)*cos(Psi)) - s*cos(Theta)), m*q*(-dTheta**2*s*cos(Psi) - g*sin(Psi)*cos(Theta) - g*sin(Theta)*cos(Psi))])
	m   = array([
		[A*(n + s)**2/12 + A*(4*h**2 - 4*h*(n - s)*sin(Theta) + (n - s)**2)/4 + M*n**2 + m*q**2 - 2*m*q*s*sin(Psi) + m*s**2, m*q*(q - s*sin(Psi))],
		[m*q*(q - s*sin(Psi)), m*q**2],
		])
	n   = array([0, 0])
	if h + q*(sin(Psi)*sin(Theta) - cos(Psi)*cos(Theta)) - s*sin(Theta) < 0:
		n = array([g*m*(q*(sin(Psi)*cos(Theta) + sin(Theta)*cos(Psi))*((sin(Psi)*sin(Theta) - cos(Psi)*cos(Theta))**2 + (sin(Psi)*cos(Theta) + sin(Theta)*cos(Psi))**2) - s*cos(Theta)), g*m*q*(sin(Psi)*cos(Theta) + sin(Theta)*cos(Psi))*((sin(Psi)*sin(Theta) - cos(Psi)*cos(Theta))**2 + (sin(Psi)*cos(Theta) + sin(Theta)*cos(Psi))**2)])
	f = f + n
	return concatenate((solve(m,f),[dTheta, dPsi]))

def DifferentialEquation(t, state, param):
	[h, o, p, q, r, n, m, M, g] = param
	[thetaDot, phiDot, theta, phi] = state
	force = array([0, 0, -o*sin(theta)*Derivative(theta, t)**2 + o*cos(theta)*Derivative(theta, t, t) - (h - y(t))*sin(theta)*Derivative(theta, t, t) - (h - y(t))*cos(theta)*Derivative(theta, t)**2 + 2*sin(theta)*Derivative(theta, t)*Derivative(y(t), t)])
	mass  = array([
		[0, 0, 0],
		[0, 0, 0],
		[0, 0, cos(theta)],
		])
	constraint  = array([0, 0])
	if -r*sin(phi) + (-p - q)*sin(theta) + y(t) < 0:
		constraint = array([-g*m*(p + q)*cos(theta), -g*m*r*cos(phi)])
	force = force + constraint
	return concatenate((solve(mass,force),[thetaDot, phiDot]))

    def test_isneginf(self):
        a = nx.array([nx.inf, -nx.inf, nx.nan, 0.0, 3.0, -3.0])
        out = nx.zeros(a.shape, bool)
        tgt = nx.array([False, True, False, False, False, False])

        res = ufl.isneginf(a)
        assert_equal(res, tgt)
        res = ufl.isneginf(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)

        a = a.astype(np.complex)
        with assert_raises(TypeError):
            ufl.isneginf(a)

    def test_fix_with_subclass(self):
        class MyArray(nx.ndarray):
            def __new__(cls, data, metadata=None):
                res = nx.array(data, copy=True).view(cls)
                res.metadata = metadata
                return res
            def __array_wrap__(self, obj, context=None):
                obj.metadata = self.metadata
                return obj

        a = nx.array([1.1, -1.1])
        m = MyArray(a, metadata='foo')
        f = ufl.fix(m)
        assert_array_equal(f, nx.array([1,-1]))
        assert_(isinstance(f, MyArray))
        assert_equal(f.metadata, 'foo')

def load(filepath, pos_filepath):
    """
    :param filepath: fits image path
    :return: an Image
    """

    # load image data information
    data = pyfits.getdata(filepath, hdu=0)
    primary = Header(pyfits.getheader(filepath, 0))
    headers = [primary]
    extcount = int(primary.get("NEXTEND", 0))

    for idx in range(1, extcount):
        ext = Header(pyfits.getheader(filepath, idx))
        headers.append(ext)

    # load position information
    pos_primary = Header(pyfits.getheader(pos_filepath, 0))
    pos_headers = [pos_primary]
    pos_extcount = int(pos_primary.get("NEXTEND", 0))

    for idx in range(1, pos_extcount):
        ext = Header(pyfits.getheader(pos_filepath, idx))
        pos_headers.append(ext)

    return Image(array(data), headers, pos_headers)

def ihfft(a, n=None, axis=-1, norm=None):
    """
    Compute the inverse FFT of a signal that has Hermitian symmetry.

    Parameters
    ----------
    a : array_like
        Input array.
    n : int, optional
        Length of the inverse FFT, the number of points along
        transformation axis in the input to use.  If `n` is smaller than
        the length of the input, the input is cropped.  If it is larger,
        the input is padded with zeros. If `n` is not given, the length of
        the input along the axis specified by `axis` is used.
    axis : int, optional
        Axis over which to compute the inverse FFT. If not given, the last
        axis is used.
    norm : {None, "ortho"}, optional
        Normalization mode (see `numpy.fft`). Default is None.

        .. versionadded:: 1.10.0

    Returns
    -------
    out : complex ndarray
        The truncated or zero-padded input, transformed along the axis
        indicated by `axis`, or the last one if `axis` is not specified.
        The length of the transformed axis is ``n//2 + 1``.

    See also
    --------
    hfft, irfft

    Notes
    -----
    `hfft`/`ihfft` are a pair analogous to `rfft`/`irfft`, but for the
    opposite case: here the signal has Hermitian symmetry in the time
    domain and is real in the frequency domain. So here it's `hfft` for
    which you must supply the length of the result if it is to be odd:

    * even: ``ihfft(hfft(a, 2*len(a) - 2) == a``, within roundoff error,
    * odd: ``ihfft(hfft(a, 2*len(a) - 1) == a``, within roundoff error.

    Examples
    --------
    >>> spectrum = np.array([ 15, -4, 0, -1, 0, -4])
    >>> np.fft.ifft(spectrum)
    array([ 1.+0.j,  2.-0.j,  3.+0.j,  4.+0.j,  3.+0.j,  2.-0.j])
    >>> np.fft.ihfft(spectrum)
    array([ 1.-0.j,  2.-0.j,  3.-0.j,  4.-0.j])

    """
    # The copy may be required for multithreading.
    a = array(a, copy=True, dtype=float)
    if n is None:
        n = a.shape[axis]
    unitary = _unitary(norm)
    output = conjugate(rfft(a, n, axis))
    return output * (1 / (sqrt(n) if unitary else n))

  def __init__(self, inc):
#    pdb.set_trace()
    self.lastUsedNumpy = 1 # Start out toggling to slow call
    self.increment = inc
    self.seed = 37
    rstate = numpy.random.get_state() # Backup the generator state.
    numpy.random.seed(37) # Set the generator state.
    alpha = array([1.0, 2.0, 1.0, 5.0, 7.0, 9.0, 3.0, 1.0])
    self.x = rdirichlet(alpha / numpy.sum(alpha))
    numpy.random.set_state(rstate) # Restore the backup state.

def schwarz_christoffel_coeff(points):
  a = exp(2j*pi*linspace(0, 1, len(points), endpoint=False))
  a.shape = (1, -1)

  p = [points[-1]] + points + points[0:1]
  b = array([ (angle(  (p[k-1]-p[k])/(p[k+1]-p[k]) )/pi)%2.0 - 1.0
              for k in xrange(1, len(p)-1) ])
  b.shape = (1,-1)

  return (a,b)