def test_resid():
    # Data is projected onto k=10 dimensional subspace then has its mean
    # removed.  Should still have rank 10.
    k = 10
    ncomp = 5
    ntotal = k
    X = np.random.standard_normal((data['nimages'], k))
    p = pca(data['fmridata'], -1, ncomp=ncomp, design_resid=X)
    yield assert_equal(
        p['basis_vectors'].shape,
        (data['nimages'], ntotal))
    yield assert_equal(
        p['basis_projections'].shape,
        data['mask'].shape + (ncomp,))
    yield assert_equal(p['pcnt_var'].shape, (ntotal,))
    yield assert_almost_equal(p['pcnt_var'].sum(), 100.)
    # if design_resid is None, we do not remove the mean, and we get
    # full rank from our data
    p = pca(data['fmridata'], -1, design_resid=None)
    rank = p['basis_vectors'].shape[1]
    yield assert_equal(rank, data['nimages'])
    rarr = reconstruct(p['basis_vectors'], p['basis_projections'], -1)
    # add back the sqrt MSE, because we standardized
    rmse = root_mse(data['fmridata'], axis=-1)[...,None]
    yield assert_array_almost_equal(rarr * rmse, data['fmridata'])

def test_scaling_io_dtype():
    # Does data dtype get set?
    # Is scaling correctly applied?
    rng = np.random.RandomState(19660520) # VBD
    ulp1_f32 = np.finfo(np.float32).eps
    types = (np.uint8, np.uint16, np.int16, np.int32, np.float32)
    with InTemporaryDirectory():
        for in_type in types:
            for out_type in types:
                data, _ = randimg_in2out(rng, in_type, out_type, 'img.nii')
                img = load_image('img.nii')
                # Check the output type is as expected
                hdr = img.metadata['header']
                assert_equal(hdr.get_data_dtype().type, out_type)
                # Check the data is within reasonable bounds. The exact bounds
                # are a little annoying to calculate - see
                # nibabel/tests/test_round_trip for inspiration
                data_back = img.get_data().copy() # copy to detach from file
                del img
                top = np.abs(data - data_back)
                nzs = (top !=0) & (data !=0)
                abs_err = top[nzs]
                if abs_err.size != 0: # all exact, that's OK.
                    continue
                rel_err = abs_err / data[nzs]
                if np.dtype(out_type).kind in 'iu':
                    slope, inter = hdr.get_slope_inter()
                    abs_err_thresh = slope / 2.0
                    rel_err_thresh = ulp1_f32
                elif np.dtype(out_type).kind == 'f':
                    abs_err_thresh = big_bad_ulp(data.astype(out_type))[nzs]
                    rel_err_thresh = ulp1_f32
                assert_true(np.all(
                    (abs_err <= abs_err_thresh) |
                    (rel_err <= rel_err_thresh)))

def test_series_from_mask():
    """ Test the smoothing of the timeseries extraction
    """
    # A delta in 3D
    data = np.zeros((40, 40, 40, 2))
    data[20, 20, 20] = 1
    mask = np.ones((40, 40, 40), dtype=np.bool)
    with InTemporaryDirectory():
        for affine in (np.eye(4), np.diag((1, 1, -1, 1)),
                        np.diag((.5, 1, .5, 1))):
            img = nib.Nifti1Image(data, affine)
            nib.save(img, 'testing.nii')
            series, header = series_from_mask('testing.nii', mask, smooth=9)
            series = np.reshape(series[:, 0], (40, 40, 40))
            vmax = series.max()
            # We are expecting a full-width at half maximum of
            # 9mm/voxel_size:
            above_half_max = series > .5*vmax
            for axis in (0, 1, 2):
                proj = np.any(np.any(np.rollaxis(above_half_max,
                                axis=axis), axis=-1), axis=-1)
                assert_equal(proj.sum(), 9/np.abs(affine[axis, axis]))

        # Check that NaNs in the data do not propagate
        data[10, 10, 10] = np.NaN
        img = nib.Nifti1Image(data, affine)
        nib.save(img, 'testing.nii')
        series, header = series_from_mask('testing.nii', mask, smooth=9)
        assert_true(np.all(np.isfinite(series)))

def test_Rintercept():
    x = F.Term('x')
    y = F.Term('x')
    xf = x.formula
    yf = y.formula
    newf = (xf+F.I)*(yf+F.I)
    assert_equal(set(newf.terms), set([x,y,x*y,sympy.Number(1)]))

def test_largest_cc():
    """ Check the extraction of the largest connected component.
    """
    a = np.zeros((6, 6, 6))
    a[1:3, 1:3, 1:3] = 1
    assert_equal(a, largest_cc(a))
    b = a.copy()
    b[5, 5, 5] = 1
    assert_equal(a, largest_cc(b))

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_clamping(I, thI=0.0, clI=256):
    regie = IconicRegistration(I, I, bins=clI)
    Ic = regie._source
    Ic2 = regie._target[1:I.shape[0]+1,1:I.shape[1]+1,1:I.shape[2]+1]
    assert_equal(Ic, Ic2.squeeze())
    dyn = Ic.max() + 1
    assert_equal(dyn, regie._joint_hist.shape[0])
    assert_equal(dyn, regie._joint_hist.shape[1])
    assert_equal(dyn, regie._source_hist.shape[0])
    assert_equal(dyn, regie._target_hist.shape[0])
    return Ic, Ic2

def test_save1():
    # A test to ensure that when a file is saved, the affine and the
    # data agree. This image comes from a NIFTI file
    img = load_image(funcfile)
    with InTemporaryDirectory():
        save_image(img, TMP_FNAME)
        img2 = load_image(TMP_FNAME)
        assert_array_almost_equal(img.affine, img2.affine)
        assert_equal(img.shape, img2.shape)
        assert_array_almost_equal(img2.get_data(), img.get_data())
        del img2

def test_update_labels():
    prng = np.random.RandomState(10)
    data, XYZ, XYZvol, vardata, signal = make_data(n=20, dim=20, r=3, 
                    amplitude=5, noise=1, jitter=1, prng=prng)
    P = os.multivariate_stat(data, vardata, XYZ)
    P.init_hidden_variables()
    p = P.data.shape[1]
    P.labels_prior = np.ones((1, p), float)
    P.label_values = np.zeros((1, p), int)
    P.labels_prior_mask = np.arange(p)
    P.update_labels()
    assert_equal(max(abs(P.labels - np.zeros(p, int))), 0)

def test_mod():
    from nipy.core.reference.coordinate_map import _matching_orth_dim
    aff = np.diag([1,2,3,1])
    for i in range(3):
        yield assert_equal(_matching_orth_dim(i, aff), (i, ''))
    aff = np.diag([-1,-2,-3,1])
    for i in range(3):
        yield assert_equal(_matching_orth_dim(i, aff), (i, ''))
    aff = np.ones((4,4))
    for i in range(3):
        val, msg = _matching_orth_dim(i, aff)
        yield assert_equal(val, None)
    aff = np.zeros((3,3))
    for i in range(2):
        val, msg = _matching_orth_dim(i, aff)
        yield assert_equal(val, None)
    aff = np.array([[1, 0, 0, 1],
                    [0, 0, 2, 1],
                    [0, 3, 0, 1],
                    [0, 0, 0, 1]])
    val, msg = _matching_orth_dim(1, aff)
    yield assert_equal(_matching_orth_dim(0, aff), (0, ''))
    yield assert_equal(_matching_orth_dim(1, aff), (2, ''))
    yield assert_equal(_matching_orth_dim(2, aff), (1, ''))
    aff = np.diag([1, 2, 0, 1])
    aff[:,3] = 1
    yield assert_equal(_matching_orth_dim(2, aff), (2, ''))

def test_save2():
    # A test to ensure that when a file is saved, the affine and the
    # data agree. This image comes from a NIFTI file 
    shape = (13,5,7,3)
    step = np.array([3.45,2.3,4.5,6.93])
    cmap = api.AffineTransform.from_start_step('ijkt', 'xyzt', [1,3,5,0], step)
    data = np.random.standard_normal(shape)
    img = api.Image(data, cmap)
    with InTemporaryDirectory():
        save_image(img, TMP_FNAME)
        img2 = load_image(TMP_FNAME)
        assert_array_almost_equal(img.affine, img2.affine)
        assert_equal(img.shape, img2.shape)
        assert_array_almost_equal(img2.get_data(), img.get_data())
        del img2

def test_names():
    # Check that the design column names are what we expect
    X = twoway.design(D, return_float=False)
    assert_equal(
        set(X.dtype.names),
        set(
            (
                "Duration_1*Weight_1",
                "Duration_1*Weight_2",
                "Duration_1*Weight_3",
                "Duration_2*Weight_1",
                "Duration_2*Weight_2",
                "Duration_2*Weight_3",
            )
        ),
    )

def test_calling_shapes():
    cs2d = CS("ij")
    cs1d = CS("i")
    cm2d = CoordinateMap(cs2d, cs2d, lambda x: x + 1)
    cm1d2d = CoordinateMap(cs1d, cs2d, lambda x: np.concatenate((x, x), axis=-1))
    at2d = AffineTransform(cs2d, cs2d, np.array([[1, 0, 1], [0, 1, 1], [0, 0, 1]]))
    at1d2d = AffineTransform(cs1d, cs2d, np.array([[1, 0], [0, 1], [0, 1]]))
    # test coordinate maps and affine transforms
    for xfm2d, xfm1d2d in ((cm2d, cm1d2d), (at2d, at1d2d)):
        arr = np.array([0, 1])
        yield assert_array_equal(xfm2d(arr), [1, 2])
        # test lists work too
        res = xfm2d([0, 1])
        yield assert_array_equal(res, [1, 2])
        # and return arrays (by checking shape attribute)
        yield assert_equal(res.shape, (2,))
        # maintaining input shape
        arr_long = arr[None, None, :]
        yield assert_array_equal(xfm2d(arr_long), arr_long + 1)
        # wrong shape array raises error
        yield assert_raises(CoordinateSystemError, xfm2d, np.zeros((3,)))
        yield assert_raises(CoordinateSystemError, xfm2d, np.zeros((3, 3)))
        # 1d to 2d
        arr = np.array(1)
        yield assert_array_equal(xfm1d2d(arr), [1, 1])
        arr_long = arr[None, None, None]
        yield assert_array_equal(xfm1d2d(arr_long), np.ones((1, 1, 2)))
        # wrong shape array raises error.  Note 1d input requires size 1
        # as final axis
        yield assert_raises(CoordinateSystemError, xfm1d2d, np.zeros((3,)))
        yield assert_raises(CoordinateSystemError, xfm1d2d, np.zeros((3, 2)))

def test_PCAMask():
    # for 2 and 4D case
    ntotal = data['nimages'] - 1
    ncomp = 5
    arr4d = data['fmridata']
    mask3d = data['mask']
    arr2d = arr4d.reshape((-1, data['nimages']))
    mask1d = mask3d.reshape((-1))
    for arr, mask in (arr4d, mask3d), (arr2d, mask1d):
        p = pca(arr, -1, mask, ncomp=ncomp)
        yield assert_equal(
            p['basis_vectors'].shape,
            (data['nimages'], ntotal))
        yield assert_equal(
            p['basis_projections'].shape,
            mask.shape + (ncomp,))
        yield assert_equal(p['pcnt_var'].shape, (ntotal,))
        yield assert_almost_equal(p['pcnt_var'].sum(), 100.)

def test_PCANoMask_nostandardize():
    ntotal = data['nimages'] - 1
    ncomp = 5
    p = pca(data['fmridata'], -1, ncomp=ncomp, standardize=False)
    assert_equal(p['basis_vectors'].shape, (data['nimages'], ntotal))
    assert_equal(p['basis_projections'].shape,  data['mask'].shape + (ncomp,))
    assert_equal(p['pcnt_var'].shape, (ntotal,))
    assert_almost_equal(p['pcnt_var'].sum(), 100.)

def test_2D():
    # check that a standard 2D PCA works too
    M = 100
    N = 20
    L = M-1 # rank after mean removal
    data = np.random.uniform(size=(M, N))
    p = pca(data)
    ts = p['basis_vectors']
    imgs = p['basis_projections']
    yield assert_equal(ts.shape, (M, L))
    yield assert_equal(imgs.shape, (L, N))
    rimgs = reconstruct(ts, imgs)
    # add back the sqrt MSE, because we standardized
    data_mean = data.mean(0)[None,...]
    demeaned = data - data_mean
    rmse = root_mse(demeaned, axis=0)[None,...]
    # also add back the mean
    yield assert_array_almost_equal((rimgs * rmse) + data_mean, data)
    # if standardize is set, or not, covariance is diagonal
    yield assert_true(diagonal_covariance(imgs))
    p = pca(data, standardize=False)
    imgs = p['basis_projections']
    yield assert_true(diagonal_covariance(imgs))

def test_input_effects():
    ntotal = data['nimages'] - 1
    # return full rank - mean PCA over last axis
    p = pca(data['fmridata'], -1)
    yield assert_equal(
        p['basis_vectors'].shape,
        (data['nimages'], ntotal))
    yield assert_equal(
        p['basis_projections'].shape,
        data['mask'].shape + (ntotal,))
    yield assert_equal(p['pcnt_var'].shape, (ntotal,))
    # Reconstructed data lacks only mean
    rarr = reconstruct(p['basis_vectors'], p['basis_projections'], -1)
    rarr = rarr + data['fmridata'].mean(-1)[...,None]
    # same effect if over axis 0, which is the default
    arr = data['fmridata']
    arr = np.rollaxis(arr, -1)
    pr = pca(arr)
    out_arr = np.rollaxis(pr['basis_projections'], 0, 4)
    yield assert_array_equal(out_arr, p['basis_projections'])
    yield assert_array_equal(p['basis_vectors'], pr['basis_vectors'])
    yield assert_array_equal(p['pcnt_var'], pr['pcnt_var'])
    # Check axis None raises error
    yield assert_raises(ValueError, pca, data['fmridata'], None)

def test_save2b():
    # A test to ensure that when a file is saved, the affine and the
    # data agree. This image comes from a NIFTI file This example has
    # a non-diagonal affine matrix for the spatial part, but is
    # 'diagonal' for the space part.  this should raise a warnings
    # about 'non-diagonal' affine matrix

    # make a 5x5 transformation
    step = np.array([3.45,2.3,4.5,6.9])
    A = np.random.standard_normal((4,4))
    B = np.diag(list(step)+[1])
    B[:4,:4] = A
    shape = (13,5,7,3)
    cmap = api.AffineTransform.from_params('ijkt', 'xyzt', B)
    data = np.random.standard_normal(shape)
    img = api.Image(data, cmap)
    with InTemporaryDirectory():
        save_image(img, TMP_FNAME)
        img2 = load_image(TMP_FNAME)
        assert_false(np.allclose(img.affine, img2.affine))
        assert_array_almost_equal(img.affine[:3,:3], img2.affine[:3,:3])
        assert_equal(img.shape, img2.shape)
        assert_array_almost_equal(img2.get_data(), img.get_data())
        del img2

def test_image_list():
    img = load_image(funcfile)
    exp_shape = (17, 21, 3, 20)
    imglst = ImageList.from_image(img, axis=-1)
    
    # Test empty ImageList
    emplst = ImageList()
    yield assert_equal(len(emplst.list), 0)

    # Test non-image construction
    a = np.arange(10)
    yield assert_raises(ValueError, ImageList, a)
    yield assert_raises(ValueError, ImageList.from_image, img, None)

    # check all the axes
    for i in range(4):
        order = range(4)
        order.remove(i)
        order.insert(0,i)
        img_re_i = img.reordered_reference(order).reordered_axes(order)
        imglst_i = ImageList.from_image(img, axis=i)

        yield assert_equal(imglst_i.list[0].shape, img_re_i.shape[1:])
        
        # check the affine as well

        yield assert_almost_equal(imglst_i.list[0].affine, 
                                  img_re_i.affine[1:,1:])

    yield assert_equal(img.shape, exp_shape)

    # length of image list should match number of frames
    yield assert_equal(len(imglst.list), img.shape[3])

    # check the affine
    A = np.identity(4)
    A[:3,:3] = img.affine[:3,:3]
    A[:3,-1] = img.affine[:3,-1]
    yield assert_almost_equal(imglst.list[0].affine, A)

    # Slicing an ImageList should return an ImageList
    sublist = imglst[2:5]
    yield assert_true(isinstance(sublist, ImageList))
    # Except when we're indexing one element
    yield assert_true(isinstance(imglst[0], Image))
    # Verify array interface
    # test __array__
    yield assert_true(isinstance(np.asarray(sublist), np.ndarray))
    # Test __setitem__
    sublist[2] = sublist[0]
    yield assert_equal(np.asarray(sublist[0]).mean(),
                       np.asarray(sublist[2]).mean())
    # Test iterator
    for x in sublist:
        yield assert_true(isinstance(x, Image))
        yield assert_equal(x.shape, exp_shape[:3])

def _test_clamping(I, thI=0.0, clI=256):
    R = HistogramRegistration(I, I, from_bins=clI)
    R.subsample(spacing=[1,1,1])
    Ic = R._from_data
    Ic2 = R._to_data[1:-1,1:-1,1:-1]
    assert_equal(Ic, Ic2)
    dyn = Ic.max() + 1
    assert_equal(dyn, R._joint_hist.shape[0])
    assert_equal(dyn, R._joint_hist.shape[1])
    return Ic, Ic2