def test_equivalent():
    ijk = CoordinateSystem("ijk")
    xyz = CoordinateSystem("xyz")
    T = np.random.standard_normal((4, 4))
    T[-1] = [0, 0, 0, 1]
    A = AffineTransform(ijk, xyz, T)

    # now, cycle through
    # all possible permutations of
    # 'ijk' and 'xyz' and confirm that
    # the mapping is equivalent

    yield assert_false, equivalent(A, A.renamed_domain({"i": "foo"}))

    try:
        import itertools

        for pijk in itertools.permutations("ijk"):
            for pxyz in itertools.permutations("xyz"):
                B = A.reordered_domain(pijk).reordered_range(pxyz)
                yield assert_true, equivalent(A, B)
    except ImportError:
        # just do some if we can't find itertools
        for pijk in ["ikj", "kij"]:
            for pxyz in ["xzy", "yxz"]:
                B = A.reordered_domain(pijk).reordered_range(pxyz)
                yield assert_true, equivalent(A, B)

def test_drop_io_dim():
    # test ordinary case of 4d to 3d
    cm4d = AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,2,3,4,1]))
    cm3d = drop_io_dim(cm4d, 't')
    yield assert_array_equal(cm3d.affine, np.diag([1, 2, 3, 1]))
    # 3d to 2d
    cm3d = AffineTransform.from_params('ijk', 'xyz', np.diag([1,2,3,1]))
    cm2d = drop_io_dim(cm3d, 'z')
    yield assert_array_equal(cm2d.affine, np.diag([1, 2, 1]))
    # test zero scaling for dropped dimension
    cm3d = AffineTransform.from_params('ijk', 'xyz', np.diag([1, 2, 0, 1]))
    cm2d = drop_io_dim(cm3d, 'z')
    yield assert_array_equal(cm2d.affine, np.diag([1, 2, 1]))
    # test not diagonal but orthogonal
    aff = np.array([[1, 0, 0, 0],
                    [0, 0, 2, 0],
                    [0, 3, 0, 0],
                    [0, 0, 0, 1]])
    cm3d = AffineTransform.from_params('ijk', 'xyz', aff)
    cm2d = drop_io_dim(cm3d, 'z')
    yield assert_array_equal(cm2d.affine, np.diag([1, 2, 1]))
    cm2d = drop_io_dim(cm3d, 'k')
    yield assert_array_equal(cm2d.affine, np.diag([1, 3, 1]))
    # and with zeros scaling for orthogonal dropped dimension
    aff[2] = 0
    cm3d = AffineTransform.from_params('ijk', 'xyz', aff)
    cm2d = drop_io_dim(cm3d, 'z')
    yield assert_array_equal(cm2d.affine, np.diag([1, 2, 1]))

def test_compose():
    value = np.array([[1.0, 2.0, 3.0]]).T
    aa = compose(E.a, E.a)
    yield assert_true, aa.inverse() is None
    yield assert_true, np.allclose(aa(value), 4 * value)
    ab = compose(E.a, E.b)
    yield assert_true, ab.inverse() is None
    assert_true, np.allclose(ab(value), 4 * value)
    ac = compose(E.a, E.c)
    yield assert_true, ac.inverse() is None
    yield assert_true, np.allclose(ac(value), value)
    bb = compose(E.b, E.b)
    #    yield assert_true, bb.inverse() is not None
    aff1 = np.diag([1, 2, 3, 1])
    affine1 = AffineTransform.from_params("ijk", "xyz", aff1)
    aff2 = np.diag([4, 5, 6, 1])
    affine2 = AffineTransform.from_params("xyz", "abc", aff2)
    # compose mapping from 'ijk' to 'abc'
    compcm = compose(affine2, affine1)
    yield assert_equal, compcm.function_domain.coord_names, ("i", "j", "k")
    yield assert_equal, compcm.function_range.coord_names, ("a", "b", "c")
    yield assert_equal, compcm.affine, np.dot(aff2, aff1)
    # check invalid coordinate mappings
    yield assert_raises, ValueError, compose, affine1, affine2

    yield assert_raises, ValueError, compose, affine1, "foo"

    cm1 = CoordinateMap(CoordinateSystem("ijk"), CoordinateSystem("xyz"), np.log)
    cm2 = CoordinateMap(CoordinateSystem("xyz"), CoordinateSystem("abc"), np.exp)
    yield assert_raises, ValueError, compose, cm1, cm2

def test_equivalent():
    ijk = CoordinateSystem('ijk')
    xyz = CoordinateSystem('xyz')
    T = np.random.standard_normal((4,4))
    T[-1] = [0,0,0,1]
    A = AffineTransform(ijk, xyz, T)

    # now, cycle through
    # all possible permutations of
    # 'ijk' and 'xyz' and confirm that
    # the mapping is equivalent

    yield assert_false, equivalent(A, A.renamed_domain({'i':'foo'}))

    try:
        import itertools
        for pijk in itertools.permutations('ijk'):
            for pxyz in itertools.permutations('xyz'):
                B = A.reordered_domain(pijk).reordered_range(pxyz)
                yield assert_true, equivalent(A, B)
    except (ImportError, AttributeError):
        # just do some if we can't find itertools, or if itertools
        # doesn't have permutations
        for pijk in ['ikj', 'kij']:
            for pxyz in ['xzy', 'yxz']:
                B = A.reordered_domain(pijk).reordered_range(pxyz)
                yield assert_true, equivalent(A, B)

def test_synchronized_order():

    data = np.random.standard_normal((3,4,7,5))
    im = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,2,3,4,1])))

    im_scrambled = im.reordered_axes('iljk').reordered_reference('xtyz')
    im_unscrambled = image.synchronized_order(im_scrambled, im)
    
    yield assert_equal, im_unscrambled.coordmap, im.coordmap
    yield assert_almost_equal, im_unscrambled.get_data(), im.get_data()
    yield assert_equal, im_unscrambled, im
    yield assert_true, im_unscrambled == im
    yield assert_false, im_unscrambled != im

    # the images don't have to be the same shape

    data2 = np.random.standard_normal((3,11,9,4))
    im2 = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,2,3,4,1])))

    im_scrambled2 = im2.reordered_axes('iljk').reordered_reference('xtyz')
    im_unscrambled2 = image.synchronized_order(im_scrambled2, im)

    yield assert_equal, im_unscrambled2.coordmap, im.coordmap

    # or the same coordmap

    data3 = np.random.standard_normal((3,11,9,4))
    im3 = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,9,3,-2,1])))

    im_scrambled3 = im3.reordered_axes('iljk').reordered_reference('xtyz')
    im_unscrambled3 = image.synchronized_order(im_scrambled3, im)

    yield assert_equal, im_unscrambled3.axes, im.axes
    yield assert_equal, im_unscrambled3.reference, im.reference

def test_compose():
    value = np.array([[1., 2., 3.]]).T
    aa = compose(E.a, E.a)
    yield assert_true, aa.inverse() is None
    yield assert_true, np.allclose(aa(value), 4*value)
    ab = compose(E.a,E.b)
    yield assert_true, ab.inverse() is None
    assert_true, np.allclose(ab(value), 4*value)
    ac = compose(E.a,E.c)
    yield assert_true, ac.inverse() is None
    yield assert_true, np.allclose(ac(value), value)
    bb = compose(E.b,E.b)
    #    yield assert_true, bb.inverse() is not None
    aff1 = np.diag([1,2,3,1])
    affine1 = AffineTransform.from_params('ijk', 'xyz', aff1)
    aff2 = np.diag([4,5,6,1])
    affine2 = AffineTransform.from_params('xyz', 'abc', aff2)
    # compose mapping from 'ijk' to 'abc'
    compcm = compose(affine2, affine1)
    yield assert_equal, compcm.function_domain.coord_names, ('i', 'j', 'k')
    yield assert_equal, compcm.function_range.coord_names, ('a', 'b', 'c')
    yield assert_equal, compcm.affine, np.dot(aff2, aff1)
    # check invalid coordinate mappings
    yield assert_raises, ValueError, compose, affine1, affine2

    yield assert_raises, ValueError, compose, affine1, 'foo'

    cm1 = CoordinateMap(CoordinateSystem('ijk'),
                        CoordinateSystem('xyz'), np.log)
    cm2 = CoordinateMap(CoordinateSystem('xyz'),
                        CoordinateSystem('abc'), np.exp)
    yield assert_raises, ValueError, compose, cm1, cm2

def test_affine_inverse():
    incs, outcs, aff = affine_v2w()
    inv = np.linalg.inv(aff)
    cm = AffineTransform(incs, outcs, aff)
    x = np.array([10, 20, 30], np.float)
    x_roundtrip = cm(cm.inverse()(x))
    yield assert_equal, x_roundtrip, x
    badaff = np.array([[1, 2, 3], [0, 0, 1]])
    badcm = AffineTransform(CoordinateSystem("ij"), CoordinateSystem("x"), badaff)
    yield assert_equal, badcm.inverse(), None

def test_rollaxis():
    data = np.random.standard_normal((3,4,7,5))
    im = Image(data, AffineTransform.from_params('ijkl', 'xyzt', np.diag([1,2,3,4,1])))

    # for the inverse we must specify an integer
    yield assert_raises, ValueError, image.rollaxis, im, 'i', True

    # Check that rollaxis preserves diagonal affines, as claimed

    yield assert_almost_equal, image.rollaxis(im, 1).affine, np.diag([2,1,3,4,1])
    yield assert_almost_equal, image.rollaxis(im, 2).affine, np.diag([3,1,2,4,1])
    yield assert_almost_equal, image.rollaxis(im, 3).affine, np.diag([4,1,2,3,1])

    # Check that ambiguous axes raise an exception
    # 'l' appears both as an axis and a reference coord name
    # and in different places

    im_amb = Image(data, AffineTransform.from_params('ijkl', 'xylt', np.diag([1,2,3,4,1])))
    yield assert_raises, ValueError, image.rollaxis, im_amb, 'l'

    # But if it's unambiguous, then
    # 'l' can appear both as an axis and a reference coord name

    im_unamb = Image(data, AffineTransform.from_params('ijkl', 'xyzl', np.diag([1,2,3,4,1])))
    im_rolled = image.rollaxis(im_unamb, 'l')
    yield assert_almost_equal, im_rolled.get_data(), \
        im_unamb.get_data().transpose([3,0,1,2])

    for i, o, n in zip('ijkl', 'xyzt', range(4)):
        im_i = image.rollaxis(im, i)
        im_o = image.rollaxis(im, o)
        im_n = image.rollaxis(im, n)

        yield assert_almost_equal, im_i.get_data(), \
                                  im_o.get_data()

        yield assert_almost_equal, im_i.affine, \
            im_o.affine

        yield assert_almost_equal, im_n.get_data(), \
            im_o.get_data()

        for _im in [im_n, im_o, im_i]:
            im_n_inv = image.rollaxis(_im, n, inverse=True)

            yield assert_almost_equal, im_n_inv.affine, \
                im.affine

            yield assert_almost_equal, im_n_inv.get_data(), \
                im.get_data()

def test__eq__():
    yield assert_true, E.a == E.a
    yield assert_false, E.a != E.a

    yield assert_false, E.a == E.b
    yield assert_true, E.a != E.b

    yield assert_true, E.singular == E.singular
    yield assert_false, E.singular != E.singular

    A = AffineTransform.from_params('ijk', 'xyz', np.diag([4,3,2,1]))
    B = AffineTransform.from_params('ijk', 'xyz', np.diag([4,3,2,1]))

    yield assert_true, A == B
    yield assert_false, A != B

def fromarray(data, innames, outnames, coordmap=None):
    """Create an image from a numpy array.

    Parameters
    ----------
    data : numpy array
        A numpy array of three dimensions.
    innames : sequence
       a list of input axis names
    innames : sequence
       a list of output axis names
    coordmap : A `CoordinateMap`
        If not specified, an identity coordinate map is created.

    Returns
    -------
    image : An `Image` object

    See Also
    --------
    load : function for loading images
    save : function for saving images

    """
    ndim = len(data.shape)
    if not coordmap:
        coordmap = AffineTransform.from_start_step(innames,
                                                   outnames,
                                                   (0.,)*ndim,
                                                   (1.,)*ndim)
                                          
    return Image(data, coordmap)

def test_affine_start_step():
    incs, outcs, aff = affine_v2w()
    start = aff[:3, 3]
    step = aff.diagonal()[:3]
    cm = AffineTransform.from_start_step(incs.coord_names, outcs.coord_names, start, step)
    yield assert_equal, cm.affine, aff
    yield assert_raises, ValueError, AffineTransform.from_start_step, "ijk", "xy", start, step

def test_renamed():

    A = AffineTransform.from_params("ijk", "xyz", np.identity(4))

    ijk = CoordinateSystem("ijk")
    xyz = CoordinateSystem("xyz")
    C = CoordinateMap(ijk, xyz, np.log)

    for B in [A, C]:
        B_re = B.renamed_domain({"i": "foo"})
        yield assert_equal, B_re.function_domain.coord_names, ("foo", "j", "k")

        B_re = B.renamed_domain({"i": "foo", "j": "bar"})
        yield assert_equal, B_re.function_domain.coord_names, ("foo", "bar", "k")

        B_re = B.renamed_range({"y": "foo"})
        yield assert_equal, B_re.function_range.coord_names, ("x", "foo", "z")

        B_re = B.renamed_range({0: "foo", 1: "bar"})
        yield assert_equal, B_re.function_range.coord_names, ("foo", "bar", "z")

        B_re = B.renamed_domain({0: "foo", 1: "bar"})
        yield assert_equal, B_re.function_domain.coord_names, ("foo", "bar", "k")

        B_re = B.renamed_range({"y": "foo", "x": "bar"})
        yield assert_equal, B_re.function_range.coord_names, ("bar", "foo", "z")

        yield assert_raises, ValueError, B.renamed_range, {"foo": "y"}
        yield assert_raises, ValueError, B.renamed_domain, {"foo": "y"}

def test_renamed():

    A = AffineTransform.from_params('ijk', 'xyz', np.identity(4))

    ijk = CoordinateSystem('ijk')
    xyz = CoordinateSystem('xyz')
    C = CoordinateMap(ijk, xyz, np.log)

    for B in [A,C]:
        B_re = B.renamed_domain({'i':'foo'})
        yield assert_equal, B_re.function_domain.coord_names, ('foo', 'j', 'k')

        B_re = B.renamed_domain({'i':'foo','j':'bar'})
        yield assert_equal, B_re.function_domain.coord_names, ('foo', 'bar', 'k')

        B_re = B.renamed_range({'y':'foo'})
        yield assert_equal, B_re.function_range.coord_names, ('x', 'foo', 'z')

        B_re = B.renamed_range({0:'foo',1:'bar'})
        yield assert_equal, B_re.function_range.coord_names, ('foo', 'bar', 'z')

        B_re = B.renamed_domain({0:'foo',1:'bar'})
        yield assert_equal, B_re.function_domain.coord_names, ('foo', 'bar', 'k')

        B_re = B.renamed_range({'y':'foo','x':'bar'})
        yield assert_equal, B_re.function_range.coord_names, ('bar', 'foo', 'z')

        yield assert_raises, ValueError, B.renamed_range, {'foo':'y'}
        yield assert_raises, ValueError, B.renamed_domain, {'foo':'y'}

def test_affine_identity():
    aff = AffineTransform.identity('ijk')
    yield assert_equal, aff.affine, np.eye(4)
    yield assert_equal, aff.function_domain, aff.function_range
    x = np.array([3, 4, 5])
    # AffineTransform's aren't CoordinateMaps, so
    # they don't have "function" attributes
    yield assert_false, hasattr(aff, 'function')

def test_product():
    affine1 = AffineTransform.from_params("i", "x", np.diag([2, 1]))
    affine2 = AffineTransform.from_params("j", "y", np.diag([3, 1]))
    affine = product(affine1, affine2)

    cm1 = CoordinateMap(CoordinateSystem("i"), CoordinateSystem("x"), np.log)

    cm2 = CoordinateMap(CoordinateSystem("j"), CoordinateSystem("y"), np.log)
    cm = product(cm1, cm2)

    yield assert_equal, affine.function_domain.coord_names, ("i", "j")
    yield assert_equal, affine.function_range.coord_names, ("x", "y")
    yield assert_almost_equal, cm([3, 4]), np.log([3, 4])
    yield assert_almost_equal, cm.function([[3, 4], [5, 6]]), np.log([[3, 4], [5, 6]])

    yield assert_equal, affine.function_domain.coord_names, ("i", "j")
    yield assert_equal, affine.function_range.coord_names, ("x", "y")
    yield assert_equal, affine.affine, np.diag([2, 3, 1])

def test_kernel():
    # Verify the the convolution with a delta function
    # gives the correct answer.
    tol = 0.9999
    sdtol = 1.0e-8
    for x in range(6):
        shape = randint(30,60,(3,))
        ii, jj, kk = randint(11,17, (3,))

        coordmap = AffineTransform.from_start_step('ijk', 'xyz', 
                                          randint(5,20,(3,))*0.25,
                                          randint(5,10,(3,))*0.5)

        signal = np.zeros(shape)
        signal[ii,jj,kk] = 1.
        signal = Image(signal, coordmap=coordmap)
    
        kernel = LinearFilter(coordmap, shape, 
                              fwhm=randint(50,100)/10.)
        ssignal = kernel.smooth(signal)
        ssignal = np.asarray(ssignal)
        ssignal[:] *= kernel.norms[kernel.normalization]

        I = np.indices(ssignal.shape)
        I.shape = (kernel.coordmap.ndim[0], np.product(shape))
        i, j, k = I[:,np.argmax(ssignal[:].flat)]

        yield assert_equal, (i,j,k), (ii,jj,kk)

        Z = kernel.coordmap(I) - kernel.coordmap([i,j,k])

        _k = kernel(Z)
        _k.shape = ssignal.shape

        yield assert_true, (np.corrcoef(_k[:].flat, ssignal[:].flat)[0,1] > tol)

        yield assert_true, ((_k[:] - ssignal[:]).std() < sdtol)
            
        def _indices(i,j,k,axis):
            I = np.zeros((3,20))
            I[0] += i
            I[1] += j
            I[2] += k
            I[axis] += np.arange(-10,10)
            return I

        vx = ssignal[i,j,(k-10):(k+10)]
        vvx = coordmap(_indices(i,j,k,2)) - coordmap([[i],[j],[k]])
        yield assert_true, (np.corrcoef(vx, kernel(vvx))[0,1] > tol)

        vy = ssignal[i,(j-10):(j+10),k]
        vvy = coordmap(_indices(i,j,k,1)) - coordmap([[i],[j],[k]])
        yield assert_true, (np.corrcoef(vy, kernel(vvy))[0,1] > tol)

        vz = ssignal[(i-10):(i+10),j,k]
        vvz = coordmap(_indices(i,j,k,0)) - coordmap([[i],[j],[k]])
        yield assert_true, (np.corrcoef(vz, kernel(vvz))[0,1] > tol)

def setup():
    def f(x):
        return 2 * x

    def g(x):
        return x / 2.0

    x = CoordinateSystem("x", "x")
    E.a = CoordinateMap(x, x, f)
    E.b = CoordinateMap(x, x, f, inverse_function=g)
    E.c = CoordinateMap(x, x, g)
    E.d = CoordinateMap(x, x, g, inverse_function=f)
    E.e = AffineTransform.identity("ijk")
    A = np.identity(4)
    A[0:3] = np.random.standard_normal((3, 4))
    E.mapping = AffineTransform.from_params("ijk", "xyz", A)
    E.singular = AffineTransform.from_params(
        "ijk", "xyzt", np.array([[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11], [8, 9, 10, 11], [0, 0, 0, 1]])
    )

def test_kernel():
    # Verify that convolution with a delta function gives the correct
    # answer.
    tol = 0.9999
    sdtol = 1.0e-8
    for x in range(6):
        shape = randint(30,60,(3,))
        # pos of delta
        ii, jj, kk = randint(11,17, (3,))
        # random affine coordmap (diagonal and translations)
        coordmap = AffineTransform.from_start_step('ijk', 'xyz', 
                                          randint(5,20,(3,))*0.25,
                                          randint(5,10,(3,))*0.5)
        # delta function in 3D array
        signal = np.zeros(shape)
        signal[ii,jj,kk] = 1.
        signal = Image(signal, coordmap=coordmap)
        # A filter with coordmap, shape matched to image
        kernel = LinearFilter(coordmap, shape, 
                              fwhm=randint(50,100)/10.)
        # smoothed normalized 3D array
        ssignal = kernel.smooth(signal).get_data()
        ssignal[:] *= kernel.norms[kernel.normalization]
        # 3 points * signal.size array
        I = np.indices(ssignal.shape)
        I.shape = (kernel.coordmap.ndims[0], np.product(shape))
        # location of maximum in smoothed array
        i, j, k = I[:, np.argmax(ssignal[:].flat)]
        # same place as we put it before smoothing?
        assert_equal((i,j,k), (ii,jj,kk))
        # get physical points position relative to position of delta
        Z = kernel.coordmap(I.T) - kernel.coordmap([i,j,k])
        _k = kernel(Z)
        _k.shape = ssignal.shape
        assert_true((np.corrcoef(_k[:].flat, ssignal[:].flat)[0,1] > tol))
        assert_true(((_k[:] - ssignal[:]).std() < sdtol))

        def _indices(i,j,k,axis):
            I = np.zeros((3,20))
            I[0] += i
            I[1] += j
            I[2] += k
            I[axis] += np.arange(-10,10)
            return I.T

        vx = ssignal[i,j,(k-10):(k+10)]
        xformed_ijk = coordmap([i, j, k])
        vvx = coordmap(_indices(i,j,k,2)) - xformed_ijk
        assert_true((np.corrcoef(vx, kernel(vvx))[0,1] > tol))
        vy = ssignal[i,(j-10):(j+10),k]
        vvy = coordmap(_indices(i,j,k,1)) - xformed_ijk
        assert_true((np.corrcoef(vy, kernel(vvy))[0,1] > tol))
        vz = ssignal[(i-10):(i+10),j,k]
        vvz = coordmap(_indices(i,j,k,0)) - xformed_ijk
        assert_true((np.corrcoef(vz, kernel(vvz))[0,1] > tol))

def test_ArrayLikeObj():
    obj = ArrayLikeObj()
    # create simple coordmap
    xform = np.eye(4)
    coordmap = AffineTransform.from_params('xyz', 'ijk', xform)
    
    # create image form array-like object and coordmap
    img = image.Image(obj, coordmap)
    yield assert_true, img.ndim == 3
    yield assert_true, img.shape == (2,3,4)
    yield assert_true, np.allclose(np.asarray(img), 1)
    yield assert_true, np.allclose(np.asarray(img), 1)
    img[:] = 4
    yield assert_true, np.allclose(np.asarray(img), 4)

def test_append_io_dim():
    aff = np.diag([1, 2, 3, 1])
    in_dims = list("ijk")
    out_dims = list("xyz")
    cm = AffineTransform.from_params(in_dims, out_dims, aff)
    cm2 = append_io_dim(cm, "l", "t")
    yield assert_array_equal(cm2.affine, np.diag([1, 2, 3, 1, 1]))
    yield assert_equal(cm2.function_range.coord_names, out_dims + ["t"])
    yield assert_equal(cm2.function_domain.coord_names, in_dims + ["l"])
    cm2 = append_io_dim(cm, "l", "t", 9, 5)
    a2 = np.diag([1, 2, 3, 5, 1])
    a2[3, 4] = 9
    yield assert_array_equal(cm2.affine, a2)
    yield assert_equal(cm2.function_range.coord_names, out_dims + ["t"])
    yield assert_equal(cm2.function_domain.coord_names, in_dims + ["l"])
    # non square case
    aff = np.array([[2, 0, 0], [0, 3, 0], [0, 0, 1], [0, 0, 1]])
    cm = AffineTransform.from_params("ij", "xyz", aff)
    cm2 = append_io_dim(cm, "q", "t", 9, 5)
    a2 = np.array([[2, 0, 0, 0], [0, 3, 0, 0], [0, 0, 0, 1], [0, 0, 5, 9], [0, 0, 0, 1]])
    yield assert_array_equal(cm2.affine, a2)
    yield assert_equal(cm2.function_range.coord_names, list("xyzt"))
    yield assert_equal(cm2.function_domain.coord_names, list("ijq"))