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)
        assert_equal(p['basis_vectors'].shape, (data['nimages'], ntotal))
        assert_equal(p['basis_projections'].shape, mask.shape + (ncomp,))
        assert_equal(p['pcnt_var'].shape, (ntotal,))
        assert_almost_equal(p['pcnt_var'].sum(), 100.)
    # Any reasonable datatype for mask
    for dt in ([np.bool_] +
               np.sctypes['int'] +
               np.sctypes['uint'] +
               np.sctypes['float']):
        p = pca(arr4d, -1, mask3d.astype(dt), ncomp=ncomp)
        assert_equal(p['basis_vectors'].shape, (data['nimages'], ntotal))
        assert_equal(p['basis_projections'].shape, mask3d.shape + (ncomp,))
        assert_equal(p['pcnt_var'].shape, (ntotal,))
        assert_almost_equal(p['pcnt_var'].sum(), 100.)
    # Mask data shape must match
    assert_raises(ValueError, pca, arr4d, -1, mask1d)

def test_apply_affine():
    XYZ = (100*(np.random.rand(10,11,12,3)-.5)).astype('int')
    T = np.eye(4)
    T[0:3,0:3] = np.random.rand(3,3)
    T[0:3,3] = 100*(np.random.rand(3)-.5)
    _XYZ = apply_affine(inverse_affine(T), apply_affine(T, XYZ))
    assert_almost_equal(_XYZ, XYZ)

def _test_similarity_measure(simi, val):
    I = Image(make_data_int16(), dummy_affine)
    J = Image(I.data.copy(), dummy_affine)
    regie = IconicRegistration(I, J)
    regie.set_source_fov(spacing=[2,1,3])
    regie.similarity = simi
    assert_almost_equal(regie.eval(np.eye(4)), val)

def test_search2():
#    Test that the search region works.

    search = rft.IntrinsicVolumes([3,4,5])
    x = np.linspace(0.1,10,100)

    stats = [rft.Gaussian(search=search)]
    ostats = [rft.Gaussian()]

    for dfn in range(5,10):
        for dfd in [40,50,np.inf]:
            stats.append(rft.FStat(dfn=dfn, dfd=dfd, search=search))
            ostats.append(rft.FStat(dfn=dfn, dfd=dfd))
            stats.append(rft.TStat(dfd=dfd, search=search))
            ostats.append(rft.TStat(dfd=dfd))
        stats.append(rft.ChiSquared(dfn=dfn, search=search))
        ostats.append(rft.ChiSquared(dfn=dfn))

    for i in range(len(stats)):
        stat = stats[i]
        ostat = ostats[i]
        v1 = stat(x)
        v2 = 0

        for j in range(search.mu.shape[0]):
            v2 += ostat.density(x, j) * search.mu[j]
        assert_almost_equal(v1, v2)

def test_evaluate_exact():
    # without mfx nor spatial relaxation
    prng = np.random.RandomState(10)
    data, XYZ, XYZvol, vardata, signal = make_data(n=20, 
            dim=np.array([20,20,20]), r=3, amplitude=5, noise=0, 
            jitter=0, prng=prng)
    p = len(signal)
    XYZvol *= 0
    XYZvol[list(XYZ)] = np.arange(p)
    P = os.multivariate_stat(data)
    P.init_hidden_variables()
    P.evaluate(nsimu=100, burnin=100, J=[XYZvol[5, 5, 5]], 
               compute_post_mean=True, verbose=verbose)
    P.log_likelihood_values = P.compute_log_region_likelihood()
    # Verify code consistency
    Q = os.multivariate_stat(data, vardata*0, XYZ, std=0, sigma=5)
    Q.init_hidden_variables()
    Q.evaluate(nsimu=100, burnin=100, J = [XYZvol[5,5,5]], 
               compute_post_mean=True, update_spatial=False, 
               verbose=verbose)
    Q.log_likelihood_values = Q.compute_log_region_likelihood()
    yield assert_almost_equal(P.mean_m.mean(),
                              Q.mean_m.mean(),
                              int(np.log10(P.nsimu))-1)
    yield assert_almost_equal(Q.log_likelihood_values.sum(),
                              P.log_likelihood_values.sum(), 0)

def test_Ragreement():
    # This code would fit the two-way ANOVA model in R

    # X = read.table('http://www-stat.stanford.edu/~jtaylo/courses/stats191/data/kidney.table', header=T)
    # names(X)
    # X$Duration = factor(X$Duration)
    # X$Weight = factor(X$Weight)
    # lm(Days~Duration*Weight, X)
    # A = anova(lm(Days~Duration*Weight, X))

    # rA = rpy.r('A')
    rA = {
        "Df": [1, 2, 2, 54],
        "F value": [7.2147239263803673, 13.120973926380339, 1.8813266871165633, np.nan],
        "Mean Sq": [209.06666666666663, 380.21666666666584, 54.51666666666663, 28.977777777777778],
        "Pr(>F)": [0.0095871255601553771, 2.2687781292164585e-05, 0.16224035152442268, np.nan],
        "Sum Sq": [209.06666666666663, 760.43333333333169, 109.03333333333326, 1564.8],
    }

    # rn = rpy.r('rownames(A)')
    rn = ["Duration", "Weight", "Duration:Weight", "Residuals"]

    pairs = [
        (rn.index("Duration"), "Duration"),
        (rn.index("Weight"), "Weight"),
        (rn.index("Duration:Weight"), "Interaction"),
    ]

    for i, j in pairs:
        assert_almost_equal(F[j], rA["F value"][i])
        assert_almost_equal(p[j], rA["Pr(>F)"][i])
        assert_almost_equal(MS[j], rA["Mean Sq"][i])
        assert_almost_equal(df[j], rA["Df"][i])
        assert_almost_equal(SS[j], rA["Sum Sq"][i])

def _test_similarity_measure(simi, val):
    I = AffineImage(make_data_int16(), dummy_affine, 'ijk')
    J = AffineImage(I.get_data().copy(), dummy_affine, 'ijk')
    R = HistogramRegistration(I, J)
    R.subsample(spacing=[2,1,3])
    R.similarity = simi
    assert_almost_equal(R.eval(Affine()), val)

def test_model_selection_exact():
    prng = np.random.RandomState(10)
    data, XYZ, XYZvol, vardata, signal = make_data(n=30, dim=20, r=3, 
                amplitude=1, noise=0, jitter=0, prng=prng)
    labels = (signal > 0).astype(int)
    P1 = os.multivariate_stat(data, labels=labels)
    P1.init_hidden_variables()
    P1.evaluate(nsimu=100, burnin=10, verbose=verbose)
    L1 = P1.compute_log_region_likelihood()
    Prior1 = P1.compute_log_prior()
    #v, m_mean, m_var = P1.v.copy(), P1.m_mean.copy(), P1.m_var.copy()
    Post1 = P1.compute_log_posterior(nsimu=1e2, burnin=1e2, verbose=verbose)
    M1 = L1 + Prior1[:-1] - Post1[:-1]
    yield assert_almost_equal(M1.mean(), 
                              P1.compute_marginal_likelihood().mean(), 0)
    P0 = os.multivariate_stat(data, labels=labels)
    P0.network *= 0
    P0.init_hidden_variables()
    P0.evaluate(nsimu=100, burnin=100, verbose=verbose)
    L0 = P0.compute_log_region_likelihood()
    Prior0 = P0.compute_log_prior()
    Post0 = P0.compute_log_posterior(nsimu=1e2, burnin=1e2, 
                                     verbose=verbose)
    M0 = L0 + Prior0[:-1] - Post0[:-1]
    yield assert_almost_equal(M0.mean(), 
                              P0.compute_marginal_likelihood().mean(), 0)
    yield assert_true(M1[1] > M0[1])
    yield assert_true(M1[0] < M0[0])

def test_alias():
    x = F.Term("x")
    f = F.aliased_function("f", lambda x: 2 * x)
    g = F.aliased_function("g", lambda x: np.sqrt(x))
    ff = F.Formula([f(x), g(x) ** 2])
    n = F.make_recarray([2, 4, 5], "x")
    yield assert_almost_equal(ff.design(n)["f(x)"], n["x"] * 2)
    yield assert_almost_equal(ff.design(n)["g(x)**2"], n["x"])

def test_alias():
    x = F.Term('x')
    f = implemented_function('f', lambda x: 2*x)
    g = implemented_function('g', lambda x: np.sqrt(x))
    ff = F.Formula([f(x), g(x)**2])
    n = F.make_recarray([2,4,5], 'x')
    yield assert_almost_equal(ff.design(n)['f(x)'], n['x']*2)
    yield assert_almost_equal(ff.design(n)['g(x)**2'], n['x'])

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_same_basis():
    arr4d = data['fmridata']
    shp = arr4d.shape
    arr2d =  arr4d.reshape((np.prod(shp[:3]), shp[3]))
    res = pos1pca(arr2d, axis=-1)
    p1b_0 = res['basis_vectors']
    for i in range(3):
        res_again = pos1pca(arr2d, axis=-1)
        assert_almost_equal(res_again['basis_vectors'], p1b_0)

def test_scaling():
    with InTemporaryDirectory():
        for dtype_type in (np.uint8, np.uint16,
                           np.int16, np.int32,
                           np.float32):
            newdata, data = uint8_to_dtype(dtype_type, 'img.nii')
            assert_almost_equal(newdata, data)
            newdata, data = float32_to_dtype(dtype_type, 'img.nii')
            assert_almost_equal(newdata, data)

def test_search3():
    # In the Gaussian case, test that search and product give same results.
    search = rft.IntrinsicVolumes([3, 4, 5, 7])
    g1 = rft.Gaussian(search=search)
    g2 = rft.Gaussian(product=search)
    x = np.linspace(0.1, 10, 100)
    y1 = g1(x)
    y2 = g2(x)
    assert_almost_equal(y1, y2)

def test_mat2vec(): 
    mat = np.eye(4)
    tmp = np.random.rand(3,3)
    U, s, Vt = np.linalg.svd(tmp)
    U /= np.linalg.det(U)
    Vt /= np.linalg.det(Vt)
    mat[0:3,0:3] = np.dot(np.dot(U, np.diag(s)), Vt)
    T = Affine(mat)
    assert_almost_equal(T.as_affine(), mat)

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_search():
    #  Test that the search region works.
    search = rft.IntrinsicVolumes([3, 4, 5])
    x = np.linspace(0.1, 10, 100)
    stat = rft.Gaussian(search=search)
    v1 = stat(x)
    v2 = (5 * x + 4 * np.sqrt(2 * np.pi)) * np.exp(-x ** 2 / 2.0) / np.power(2 * np.pi, 1.5) + 3 * scipy.stats.norm.sf(
        x
    )
    assert_almost_equal(v1, v2)

def test_search4():
    # Test that the search/product work well together
    search = rft.IntrinsicVolumes([3, 4, 5])
    product = rft.IntrinsicVolumes([1, 2])
    x = np.linspace(0.1, 10, 100)
    g1 = rft.Gaussian()
    g2 = rft.Gaussian(product=product)
    y = g2(x, search=search)
    z = g1(x, search=search * product)
    assert_almost_equal(y, z)

def test_both():
    k1 = 10
    k2 = 8
    ncomp = 5
    ntotal = k1
    X1 = np.random.standard_normal((data['nimages'], k1))
    X2 = np.random.standard_normal((data['nimages'], k2))
    p = pca(data['fmridata'], -1, ncomp=ncomp, design_resid=X2, design_keep=X1)
    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_keep():
    # 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_keep=X)
    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.)