def group_analysis_signs(design, contrast, mask, signs=None):
    """
    This function refits the EM model with a vector of signs.
    Used in the permutation tests.

    Returns the maximum of the T-statistic within mask

    Parameters
    ----------

    design: one of 'block', 'event'

    contrast: str

    mask: array-like

    signs: ndarray, optional
         Defaults to np.ones. Should have shape (*,nsubj)
         where nsubj is the number of effects combined in the group analysis.

    Returns
    -------

    minT: np.ndarray, minima of T statistic within mask, one for each
         vector of signs

    maxT: np.ndarray, maxima of T statistic within mask, one for each
         vector of signs
    
    """

    maska = np.asarray(mask).astype(np.bool)

    # Which subjects have this (contrast, design) pair?

    subjects = futil.subject_dirs(design, contrast)

    sd = np.array([np.array(load_image(pjoin(s, "sd.nii")))[:,maska]
                   for s in subjects])
    Y = np.array([np.array(load_image(pjoin(s, "effect.nii")))[:,maska]
                  for s in subjects])

    if signs is None:
        signs = np.ones((1, Y.shape[0]))

    maxT = np.empty(signs.shape[0])
    minT = np.empty(signs.shape[0])

    for i, sign in enumerate(signs):
        signY = sign[:,np.newaxis] * Y
        varest = onesample.estimate_varatio(signY, sd)
        random_var = varest['random']

        adjusted_var = sd**2 + random_var
        adjusted_sd = np.sqrt(adjusted_var)

        results = onesample.estimate_mean(Y, adjusted_sd) 
        T = results['t']
        minT[i], maxT[i] = np.nanmin(T), np.nanmax(T)
    return minT, maxT

def group_analysis(design, contrast):
    """ Compute group analysis effect, t, sd for `design` and `contrast`

    Saves to disk in 'group' analysis directory

    Parameters
    ----------
    design : {'block', 'event'}
    contrast : str
        contrast name
    """
    array = np.array # shorthand
    # Directory where output will be written
    odir = futil.ensure_dir(futil.DATADIR, 'group', design, contrast)

    # Which subjects have this (contrast, design) pair?
    subj_con_dirs = futil.subj_des_con_dirs(design, contrast)
    if len(subj_con_dirs) == 0:
        raise ValueError('No subjects for %s, %s' % (design, contrast))

    # Assemble effects and sds into 4D arrays
    sds = []
    Ys = []
    for s in subj_con_dirs:
        sd_img = load_image(pjoin(s, "sd.nii"))
        effect_img = load_image(pjoin(s, "effect.nii"))
        sds.append(sd_img.get_data())
        Ys.append(effect_img.get_data())
    sd = array(sds)
    Y = array(Ys)

    # This function estimates the ratio of the fixed effects variance
    # (sum(1/sd**2, 0)) to the estimated random effects variance
    # (sum(1/(sd+rvar)**2, 0)) where rvar is the random effects variance.

    # The EM algorithm used is described in:
    #
    # Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H.,
    #    Morales, F., Evans, A.C. (2002). \'A general statistical
    #    analysis for fMRI data\'. NeuroImage, 15:1-15
    varest = onesample.estimate_varatio(Y, sd)
    random_var = varest['random']

    # XXX - if we have a smoother, use
    # random_var = varest['fixed'] * smooth(varest['ratio'])

    # Having estimated the random effects variance (and possibly smoothed it),
    # the corresponding estimate of the effect and its variance is computed and
    # saved.

    # This is the coordmap we will use
    coordmap = futil.load_image_fiac("fiac_00","wanatomical.nii").coordmap

    adjusted_var = sd**2 + random_var
    adjusted_sd = np.sqrt(adjusted_var)

    results = onesample.estimate_mean(Y, adjusted_sd) 
    for n in ['effect', 'sd', 't']:
        im = api.Image(results[n], copy(coordmap))
        save_image(im, pjoin(odir, "%s.nii" % n))

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)
    save_image(img, tmpfile.name)
    img2 = load_image(tmpfile.name)
    yield assert_true, np.allclose(img.affine, img2.affine)
    yield assert_equal, img.shape, img2.shape
    yield assert_true, np.allclose(np.asarray(img2), np.asarray(img))

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 setUp(self):
        self.fd = np.asarray(load_image(funcfile))
        self.fi = FmriImageList.from_image(load_image(funcfile))
        # I think it makes more sense to use fd instead of fi for GLM
        # purposes -- reduces some noticeable overhead in creating the
        # array from FmriImageList

        # create a design matrix, model and contrast matrix

        self.design = noise((self.fd.shape[0],3))
        self.model = ols_model(self.design)
        self.cmatrix = np.array([[1,0,0],[0,1,0]])

def fixed_effects(subj, design):
    """ Fixed effects (within subject) for FIAC model

    Finds run by run estimated model results, creates fixed effects results
    image per subject.

    Parameters
    ----------
    subj : int
        subject number 1..6 inclusive
    design : {'standard'}
        design type
    """
    # First, find all the effect and standard deviation images
    # for the subject and this design type
    path_dict = futil.path_info_design(subj, design)
    rootdir = path_dict['rootdir']
    # The output directory
    fixdir = pjoin(rootdir, "fixed")
    # Fetch results images from run estimations
    results = futil.results_table(path_dict)
    # Get our hands on the relevant coordmap to save our results
    coordmap = futil.load_image_fiac("_%02d" % subj,
                                     "wanatomical.nii").coordmap
    # Compute the "fixed" effects for each type of contrast
    for con in results:
        fixed_effect = 0
        fixed_var = 0
        for effect, sd in results[con]:
            effect = load_image(effect).get_data()
            sd = load_image(sd).get_data()
            var = sd ** 2

            # The optimal, in terms of minimum variance, combination of the
            # effects has weights 1 / var
            #
            # XXX regions with 0 variance are set to 0
            # XXX do we want this or np.nan?
            ivar = np.nan_to_num(1. / var)
            fixed_effect += effect * ivar
            fixed_var += ivar

        # Now, compute the fixed effects variance and t statistic
        fixed_sd = np.sqrt(fixed_var)
        isd = np.nan_to_num(1. / fixed_sd)
        fixed_t = fixed_effect * isd

        # Save the results
        odir = futil.ensure_dir(fixdir, con)
        for a, n in zip([fixed_effect, fixed_sd, fixed_t],
                        ['effect', 'sd', 't']):
            im = api.Image(a, copy(coordmap))
            save_image(im, pjoin(odir, '%s.nii' % n))

def get_fmri_anat(path_dict):
    """Get the images for a given subject/run.

    Returns
    -------
    fmri : ndarray
    
    anat : NIPY image
    """
    fmri = np.array(load_image(pjoin("%(rootdir)s/swafunctional_%(run)02d.nii") % path_dict))
    anat = load_image(pjoin(DATADIR, "fiac_%(subj)02d", "wanatomical.nii") % path_dict)
    return fmri, anat

def test_write():
    fname = "myfile.nii"
    img = load_image(funcfile)
    with InTemporaryDirectory():
        save_image(img, fname)
        test = FmriImageList.from_image(load_image(fname))
        assert_equal(test[0].affine.shape, (4, 4))
        assert_equal(img[0].affine.shape, (5, 4))
        # Check the affine...
        A = np.identity(4)
        A[:3, :3] = img[:, :, :, 0].affine[:3, :3]
        A[:3, -1] = img[:, :, :, 0].affine[:3, -1]
        assert_true(np.allclose(test[0].affine, A))
        del test

def test_space_time_realign():
    path, fname = psplit(funcfile)
    original_affine = load_image(funcfile).affine
    path, fname = psplit(funcfile)
    froot, _ = fname.split('.', 1)
    with InTemporaryDirectory():
        # Make another image with .nii extension and extra dot in filename
        save_image(load_image(funcfile), 'my.test.nii')
        for in_fname, out_fname in ((funcfile, froot + '_mc.nii.gz'),
                                    ('my.test.nii', 'my.test_mc.nii.gz')):
            xforms = reg.space_time_realign(in_fname, 2.0, out_name='.')
            assert_true(np.allclose(xforms[0].as_affine(), np.eye(4), atol=1e-7))
            assert_false(np.allclose(xforms[-1].as_affine(), np.eye(4), atol=1e-3))
            img = load_image(out_fname)
            npt.assert_almost_equal(original_affine, img.affine)

def fixed_effects(subj, design):
    """
    Fixed effects (within subject) for FIAC model
    """

    # First, find all the effect and standard deviation images
    # for the subject and this design type

    path_dict = futil.path_info2(subj, design)
    rootdir = path_dict['rootdir']
    # The output directory
    fixdir = pjoin(rootdir, "fixed")

    results = futil.results_table(path_dict)

    # Get our hands on the relevant coordmap to
    # save our results
    coordmap = futil.load_image_fiac("fiac_%02d" % subj,
                                     "wanatomical.nii").coordmap

    # Compute the "fixed" effects for each type of contrast
    for con in results:
        fixed_effect = 0
        fixed_var = 0
        for effect, sd in results[con]:
            effect = load_image(effect); sd = load_image(sd)
            var = np.array(sd)**2

            # The optimal, in terms of minimum variance, combination of the
            # effects has weights 1 / var
            #
            # XXX regions with 0 variance are set to 0
            # XXX do we want this or np.nan?
            ivar = np.nan_to_num(1. / var)
            fixed_effect += effect * ivar
            fixed_var += ivar

        # Now, compute the fixed effects variance and t statistic
        fixed_sd = np.sqrt(fixed_var)
        isd = np.nan_to_num(1. / fixed_sd)
        fixed_t = fixed_effect * isd

        # Save the results
        odir = futil.ensure_dir(fixdir, con)
        for a, n in zip([fixed_effect, fixed_sd, fixed_t],
                        ['effect', 'sd', 't']):
            im = api.Image(a, coordmap.copy())
            save_image(im, pjoin(odir, '%s.nii' % n))

def group_analysis(design, contrast):
    """
    Compute group analysis effect, sd and t
    for a given contrast and design type
    """
    array = np.array # shorthand
    # Directory where output will be written
    odir = futil.ensure_dir(futil.DATADIR, 'group', design, contrast)

    # Which subjects have this (contrast, design) pair?
    subjects = futil.subject_dirs(design, contrast)

    sd = array([array(load_image(pjoin(s, "sd.nii"))) for s in subjects])
    Y = array([array(load_image(pjoin(s, "effect.nii"))) for s in subjects])

    # This function estimates the ratio of the
    # fixed effects variance (sum(1/sd**2, 0))
    # to the estimated random effects variance
    # (sum(1/(sd+rvar)**2, 0)) where
    # rvar is the random effects variance.

    # The EM algorithm used is described in 
    #
    # Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H., 
    #    Morales, F., Evans, A.C. (2002). \'A general statistical 
    #    analysis for fMRI data\'. NeuroImage, 15:1-15

    varest = onesample.estimate_varatio(Y, sd)
    random_var = varest['random']

    # XXX - if we have a smoother, use
    # random_var = varest['fixed'] * smooth(varest['ratio'])

    # Having estimated the random effects variance (and
    # possibly smoothed it), the corresponding
    # estimate of the effect and its variance is
    # computed and saved.

    # This is the coordmap we will use
    coordmap = futil.load_image_fiac("fiac_00","wanatomical.nii").coordmap

    adjusted_var = sd**2 + random_var
    adjusted_sd = np.sqrt(adjusted_var)

    results = onesample.estimate_mean(Y, adjusted_sd) 
    for n in ['effect', 'sd', 't']:
        im = api.Image(results[n], coordmap.copy())
        save_image(im, pjoin(odir, "%s.nii" % n))

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 transformatio
    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)

    save_image(img, tmpfile.name)
    img2 = load_image(tmpfile.name)
    yield assert_false, np.allclose(img.affine, img2.affine)
    yield assert_true, np.allclose(img.affine[:3,:3], img2.affine[:3,:3])
    yield assert_equal, img.shape, img2.shape
    yield assert_true, np.allclose(np.asarray(img2), np.asarray(img))

def test_save3():
    # A test to ensure that when a file is saved, the affine
    # and the data agree. In this case, things don't agree:
    # i) the pixdim is off
    # ii) makes the affine off
    step = np.array([3.45,2.3,4.5,6.9])
    shape = (13,5,7,3)
    mni_xyz = mni_csm(3).coord_names
    cmap = AT(CS('jkli'),
              CS(('t',) + mni_xyz[::-1]),
              from_matvec(np.diag([0,3,5,1]), step))
    data = np.random.standard_normal(shape)
    img = api.Image(data, cmap)
    # with InTemporaryDirectory():
    with InTemporaryDirectory():
        save_image(img, TMP_FNAME)
        tmp = load_image(TMP_FNAME)
        # Detach image from file so we can delete it
        data = tmp.get_data().copy()
        img2 = api.Image(data, tmp.coordmap, tmp.metadata)
        del tmp
    assert_equal(tuple([img.shape[l] for l in [3,2,1,0]]), img2.shape)
    a = np.transpose(np.asarray(img), [3,2,1,0])
    assert_false(np.allclose(img.affine, img2.affine))
    assert_true(np.allclose(a, img2.get_data()))

def setup():
    tmp_img = load_image(funcfile)
    # For now, img is an Image
    # instead of an XYZImage
    A = np.identity(4)
    A[:3,:3] = tmp_img.affine[:3,:3]
    A[:3,-1] = tmp_img.affine[:3,-1]
    xyz_data = tmp_img.get_data()
    xyz_img = XYZImage(xyz_data, A,
                       tmp_img.axes.coord_names[:3] + ('t',))
    # If load_image returns an XYZImage, I'd really be
    # starting from xyz_img, so from here

    # Here, I'm just doing this so I know that
    # img.shape[0] is the number of volumes
    img = image_rollaxis(Image(xyz_img._data, xyz_img.coordmap), 't')
    data_dict['nimages'] = img.shape[0]
    # It might be worth to make
    # data a public attribute of Image/XYZImage
    # and we might rename get_data-> data_as_array = np.asarray(self.data)
    # Then, the above would not access a private attribute
    
    # Below, I am just making a mask
    # because I already have img, I 
    # know I can do this
    # In principle, though, the pca function
    # will just take another XYZImage as a mask
    img_data = img.get_data()
    first_frame = img_data[0]
    mask = XYZImage(np.greater(np.asarray(first_frame), 500).astype(np.float64), A, xyz_img.axes.coord_names[:3])
    data_dict['fmridata'] = xyz_img
    data_dict['mask'] = mask
    data_dict['img'] = img
    print data_dict['mask'].shape, np.sum(np.array(data_dict['mask']))

def permutation_test(design, contrast, mask=GROUP_MASK, nsample=1000):
    """
    Perform a permutation (sign) test for a given design type and
    contrast. It is a Monte Carlo test because we only sample nsample
    possible sign arrays.

    Parameters
    ----------
    design: one of ['block', 'event']
    contrast: str
    nsample: int

    Returns
    -------
    min_vals: np.ndarray
    max_vals: np.ndarray
    """
    maska = np.asarray(mask).astype(np.bool)
    subjects = futil.subject_dirs(design, contrast)
    Y = np.array([np.array(load_image(pjoin(s, "effect.nii")))[:,maska]
                  for s in subjects])
    nsubj = Y.shape[0]
    signs = 2*np.greater(np.random.sample(size=(nsample, nsubj)), 0.5) - 1
    min_vals, max_vals = group_analysis_signs(design, contrast, maska, signs)
    return min_vals, max_vals

def setup():
    img = load_image(funcfile)
    arr = np.array(img)
    #arr = np.rollaxis(arr, 3)
    data['nimages'] = arr.shape[3]
    data['fmridata'] = arr
    frame = data['fmridata'][...,0]
    data['mask'] = (frame > 500).astype(np.float64)

def test_labels1():
    img = load_image(funcfile)
    parcelmap = fromarray(np.asarray(img[0]), 'kji', 'zyx')    
    parcelmap = (np.asarray(parcelmap) * 100).astype(np.int32)
    v = 0
    for i, d in fmri_generator(img, parcels(parcelmap)):
        v += d.shape[1]
    assert_equal(v, parcelmap.size)

def test_model_out_img():
    # Model output image
    cmap = load_image(anatfile).coordmap
    shape = (2,3,4)
    fname = 'myfile.nii'
    with InTemporaryDirectory():
        moi = ModelOutputImage(fname, cmap, shape)
        for i in range(shape[0]):
            moi[i] = i
        for i in range(shape[0]):
            assert_array_equal(moi[i], i)
        moi.save()
        assert_raises(ValueError, moi.__setitem__, 0, 1)
        assert_raises(ValueError, moi.__getitem__, 0)
        new_img = load_image(fname)
        for i in range(shape[0]):
            assert_array_equal(new_img[i].get_data(), i)

def test_labels1():
    img = load_image(funcfile)
    data = img.get_data()
    parcelmap = Image(img[0].get_data(), AfT("kji", "zyx", np.eye(4)))
    parcelmap = (parcelmap.get_data() * 100).astype(np.int32)
    v = 0
    for i, d in axis0_generator(data, parcels(parcelmap)):
        v += d.shape[1]
    assert_equal(v, parcelmap.size)

def test_write():
    fp, fname = mkstemp('.nii')
    img = load_image(funcfile)
    save_image(img, fname)
    test = FmriImageList.from_image(load_image(fname))
    yield assert_equal, test[0].affine.shape, (4,4)
    yield assert_equal, img[0].affine.shape, (5,4)

    # Check the affine...
    A = np.identity(4)
    A[:3,:3] = img[:,:,:,0].affine[:3,:3]
    A[:3,-1] = img[:,:,:,0].affine[:3,-1]
    yield assert_true, np.allclose(test[0].affine, A)

    # Under windows, if you don't close before delete, you get a
    # locking error.
    os.close(fp)
    os.remove(fname)