How to use the tedana.utils.filewrite function in tedana

Default: None
    dindex : int, optional
        Index of data (2nd dimension) to use for clustering. Default: 0
    tindex : int, optional
        Index of data (2nd dimension) to use for thresholding. Default: 0

    Returns
    -------
    clustered : (S x T) np.ndarray
        Boolean array indicated data samples to be retained after clustering
    """

    if infile is None:
        data = data.copy()
        data[data < thr] = 0
        infile = filewrite(unmask(data, mask), '__clin.nii.gz', ref_img, gzip=True)

    addopts = ''
    if data is not None and data.squeeze().ndim > 1 and dindex + tindex == 0:
        addopts = '-doall'
    else:
        addopts = '-1dindex {0} -1tindex {1}'.format(str(dindex), str(tindex))

    cmd_str = '3dmerge -overwrite {0} -dxyz=1 -1clust 1 {1:d} ' \
              '-1thresh {2:.02f} -prefix __clout.nii.gz {3}'
    os.system(cmd_str.format(addopts, int(csize), float(thr), infile))

    clustered = load_image('__clout.nii.gz')[mask] != 0
    return clustered

lgr.info('++ Computing T2* map')
    t2s, s0, t2ss, s0s, t2sG, s0G = t2sadmap(catd, tes, mask, masksum,
                                             start_echo=1)

    # set a hard cap for the T2* map
    # anything that is 10x higher than the 99.5 %ile will be reset to 99.5 %ile
    cap_t2s = stats.scoreatpercentile(t2s.flatten(), 99.5,
                                      interpolation_method='lower')
    t2s[t2s > cap_t2s * 10] = cap_t2s
    filewrite(t2s, op.join(out_dir, 't2sv'), ref_img)
    filewrite(s0, op.join(out_dir, 's0v'), ref_img)
    filewrite(t2ss, op.join(out_dir, 't2ss'), ref_img)
    filewrite(s0s, op.join(out_dir, 's0vs'), ref_img)
    filewrite(t2sG, op.join(out_dir, 't2svG'), ref_img)
    filewrite(s0G, op.join(out_dir, 's0vG'), ref_img)

    # optimally combine data
    OCcatd = make_optcom(catd, t2sG, tes, mask, combmode)

    # regress out global signal unless explicitly not desired
    if not no_gscontrol:
        catd, OCcatd = gscontrol_raw(catd, OCcatd, n_echos, ref_img)

    if mixm is None:
        lgr.info("++ Doing ME-PCA and ME-ICA")
        n_components, dd = tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, ref_img,
                                  tes=tes, kdaw=kdaw, rdaw=rdaw, ste=ste)
        mmix_orig = tedica(n_components, dd, conv, fixed_seed, cost=initcost,
                           final_cost=finalcost)
        np.savetxt(op.join(out_dir, '__meica_mix.1D'), mmix_orig)
        seldict, comptable, betas, mmix = fitmodels_direct(catd, mmix_orig,

sphis -= sphis.mean()
    filewrite(unmask(sphis, Gmask), 'T1gs', ref_img)

    # find time course ofc the spatial global signal
    # make basis with the Legendre basis
    glsig = np.linalg.lstsq(np.atleast_2d(sphis).T, dat)[0]
    glsig = stats.zscore(glsig, axis=None)
    np.savetxt('glsig.1D', glsig)
    glbase = np.hstack([Lmix, glsig.T])

    # Project global signal out of optimally combined data
    sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T)[0]
    tsoc_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T,
                             np.atleast_2d(glbase.T[dtrank])) + Gmu[Gmask][:, np.newaxis]

    filewrite(optcom, 'tsoc_orig', ref_img)
    dm_optcom = unmask(tsoc_nogs, Gmask)
    filewrite(dm_optcom, 'tsoc_nogs', ref_img)

    # Project glbase out of each echo
    dm_catd = catd.copy()  # don't overwrite catd
    for echo in range(n_echos):
        dat = dm_catd[:, echo, :][Gmask]
        sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T)[0]
        e_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T,
                              np.atleast_2d(glbase.T[dtrank]))
        dm_catd[:, echo, :] = unmask(e_nogs, Gmask)

    return dm_catd, dm_optcom

'dn_ts_OC_T1c', ref_img)

    """
    Orthogonalize mixing matrix w.r.t. T1-GS
    """
    mmixnogs = mmix.T - np.dot(np.linalg.lstsq(glsig.T, mmix)[0].T, glsig)
    mmixnogs_mu = mmixnogs.mean(-1)
    mmixnogs_std = mmixnogs.std(-1)
    mmixnogs_norm = (mmixnogs - mmixnogs_mu[:, np.newaxis]) / mmixnogs_std[:, np.newaxis]
    mmixnogs_norm = np.vstack([np.atleast_2d(np.ones(max(glsig.shape))), glsig, mmixnogs_norm])

    """
    Write T1-GS corrected components and mixing matrix
    """
    sol = np.linalg.lstsq(mmixnogs_norm.T, dat.T)
    filewrite(unmask(sol[0].T[:, 2:], mask), 'betas_hik_OC_T1c', ref_img)
    np.savetxt('meica_mix_T1c.1D', mmixnogs)

mask, masksum = make_adaptive_mask(catd, minimum=False, getsum=True)

    lgr.info('++ Computing T2* map')
    t2s, s0, t2ss, s0s, t2sG, s0G = t2sadmap(catd, tes, mask, masksum,
                                             start_echo=1)

    # set a hard cap for the T2* map
    # anything that is 10x higher than the 99.5 %ile will be reset to 99.5 %ile
    cap_t2s = stats.scoreatpercentile(t2s.flatten(), 99.5,
                                      interpolation_method='lower')
    t2s[t2s > cap_t2s * 10] = cap_t2s
    filewrite(t2s, op.join(out_dir, 't2sv'), ref_img)
    filewrite(s0, op.join(out_dir, 's0v'), ref_img)
    filewrite(t2ss, op.join(out_dir, 't2ss'), ref_img)
    filewrite(s0s, op.join(out_dir, 's0vs'), ref_img)
    filewrite(t2sG, op.join(out_dir, 't2svG'), ref_img)
    filewrite(s0G, op.join(out_dir, 's0vG'), ref_img)

    # optimally combine data
    OCcatd = make_optcom(catd, t2sG, tes, mask, combmode)

    # regress out global signal unless explicitly not desired
    if not no_gscontrol:
        catd, OCcatd = gscontrol_raw(catd, OCcatd, n_echos, ref_img)

    if mixm is None:
        lgr.info("++ Doing ME-PCA and ME-ICA")
        n_components, dd = tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, ref_img,
                                  tes=tes, kdaw=kdaw, rdaw=rdaw, ste=ste)
        mmix_orig = tedica(n_components, dd, conv, fixed_seed, cost=initcost,
                           final_cost=finalcost)
        np.savetxt(op.join(out_dir, '__meica_mix.1D'), mmix_orig)

LGR.info('Variance explained by ICA decomposition: '
             '{:.02f}%'.format(varexpl))

    # create component and de-noised time series and save to files
    hikts = betas[:, acc].dot(mmix.T[acc, :])
    midkts = betas[:, midk].dot(mmix.T[midk, :])
    lowkts = betas[:, rej].dot(mmix.T[rej, :])
    dnts = data[mask] - lowkts - midkts

    if len(acc) != 0:
        fout = utils.filewrite(utils.unmask(hikts, mask),
                               'hik_ts_{0}'.format(suffix), ref_img)
        LGR.info('Writing high-Kappa time series: {}'.format(op.abspath(fout)))

    if len(midk) != 0:
        fout = utils.filewrite(utils.unmask(midkts, mask),
                               'midk_ts_{0}'.format(suffix), ref_img)
        LGR.info('Writing mid-Kappa time series: {}'.format(op.abspath(fout)))

    if len(rej) != 0:
        fout = utils.filewrite(utils.unmask(lowkts, mask),
                               'lowk_ts_{0}'.format(suffix), ref_img)
        LGR.info('Writing low-Kappa time series: {}'.format(op.abspath(fout)))

    fout = utils.filewrite(utils.unmask(dnts, mask),
                           'dn_ts_{0}'.format(suffix), ref_img)
    LGR.info('Writing denoised time series: {}'.format(op.abspath(fout)))

    return varexpl

# find time course ofc the spatial global signal
    # make basis with the Legendre basis
    glsig = np.linalg.lstsq(np.atleast_2d(sphis).T, dat)[0]
    glsig = stats.zscore(glsig, axis=None)
    np.savetxt('glsig.1D', glsig)
    glbase = np.hstack([Lmix, glsig.T])

    # Project global signal out of optimally combined data
    sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T)[0]
    tsoc_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T,
                             np.atleast_2d(glbase.T[dtrank])) + Gmu[Gmask][:, np.newaxis]

    filewrite(optcom, 'tsoc_orig', ref_img)
    dm_optcom = unmask(tsoc_nogs, Gmask)
    filewrite(dm_optcom, 'tsoc_nogs', ref_img)

    # Project glbase out of each echo
    dm_catd = catd.copy()  # don't overwrite catd
    for echo in range(n_echos):
        dat = dm_catd[:, echo, :][Gmask]
        sol = np.linalg.lstsq(np.atleast_2d(glbase), dat.T)[0]
        e_nogs = dat - np.dot(np.atleast_2d(sol[dtrank]).T,
                              np.atleast_2d(glbase.T[dtrank]))
        dm_catd[:, echo, :] = unmask(e_nogs, Gmask)

    return dm_catd, dm_optcom

lgr.info('++ Computing Mask')
    mask, masksum = make_adaptive_mask(catd, minimum=False, getsum=True)

    lgr.info('++ Computing T2* map')
    t2s, s0, t2ss, s0s, t2sG, s0G = t2sadmap(catd, tes, mask, masksum,
                                             start_echo=1)

    # set a hard cap for the T2* map
    # anything that is 10x higher than the 99.5 %ile will be reset to 99.5 %ile
    cap_t2s = stats.scoreatpercentile(t2s.flatten(), 99.5,
                                      interpolation_method='lower')
    t2s[t2s > cap_t2s * 10] = cap_t2s
    filewrite(t2s, op.join(out_dir, 't2sv'), ref_img)
    filewrite(s0, op.join(out_dir, 's0v'), ref_img)
    filewrite(t2ss, op.join(out_dir, 't2ss'), ref_img)
    filewrite(s0s, op.join(out_dir, 's0vs'), ref_img)
    filewrite(t2sG, op.join(out_dir, 't2svG'), ref_img)
    filewrite(s0G, op.join(out_dir, 's0vG'), ref_img)

    # optimally combine data
    OCcatd = make_optcom(catd, t2sG, tes, mask, combmode)

    # regress out global signal unless explicitly not desired
    if not no_gscontrol:
        catd, OCcatd = gscontrol_raw(catd, OCcatd, n_echos, ref_img)

    if mixm is None:
        lgr.info("++ Doing ME-PCA and ME-ICA")
        n_components, dd = tedpca(catd, OCcatd, combmode, mask, t2s, t2sG, stabilize, ref_img,
                                  tes=tes, kdaw=kdaw, rdaw=rdaw, ste=ste)
        mmix_orig = tedica(n_components, dd, conv, fixed_seed, cost=initcost,
                           final_cost=finalcost)

if full_sel:
        lgr.info('++ Performing spatial clustering of components')
        for i in range(n_components):
            # save out files
            out = np.zeros((n_samp, 4))
            if fout is not None:
                ccname, gzip = 'cc{:03d}'.format(i), False
            else:
                ccname, gzip = '.cc_temp.nii.gz', True

            out[:, 0] = np.squeeze(unmask(PSC[:, i], mask))
            out[:, 1] = np.squeeze(unmask(F_R2_maps[:, i], t2s != 0))
            out[:, 2] = np.squeeze(unmask(F_S0_maps[:, i], t2s != 0))
            out[:, 3] = np.squeeze(unmask(Z_maps[:, i], mask))

            ccname = filewrite(out, ccname, ref_img, gzip=gzip)

            if get_dtype(ref_img) == 'GIFTI':
                continue  # TODO: pass through GIFTI file data as below

            os.system('3drefit -sublabel 0 PSC -sublabel 1 F_R2 -sublabel 2 F_SO '
                      '-sublabel 3 Z_sn {} 2> /dev/null > /dev/null'.format(ccname))

            csize = np.max([int(n_voxels * 0.0005) + 5, 20])

            # Do simple clustering on F
            # TODO: can be replaced with nilearn.image.threshold_img
            # TODO: fmin is being cast to an integer here -- is that purposeful?!
            os.system('3dcalc -overwrite -a {}[1..2] -expr \'a*step(a-{:0d})\' -prefix '
                      '.fcl_in.nii.gz -overwrite'.format(ccname, int(fmin)))
            # TODO: can be replaced with nilearn.regions.connected_regions
            os.system('3dmerge -overwrite -dxyz=1 -1clust 1 {:0d} -doall '

mask : (S,) array_like
        Boolean mask array
    ref_img : str or img_like
        Reference image to dictate how outputs are saved to disk
    suffix : str, optional
        Appended to name of saved files (before extension). Default: ''

    Returns
    -------
    fname : str
        Filepath to saved file
    """

    # write feature versions of components
    feats = unmask(computefeats2(data, mmix, mask), mask)
    fname = filewrite(feats, 'feats_{0}'.format(suffix), ref_img)

    return fname