def test_sub_treatment(self):

        t = FMRITreatment(output_dir=self.tmp_dir)
        t.enable_draft_testing()
        sub_ts = t.split()
        for sub_t in sub_ts:
            sub_t.run()

 def test_default_treatment_parallel_LAN(self):
     #pyhrf.verbose.set_verbosity(1)
     if cfg['parallel-LAN']['enable_unit_test'] == 1:
         t = FMRITreatment(make_outputs=False, result_dump_file=None,
                           output_dir=self.tmp_dir)
         t.enable_draft_testing()
         t.run(parallel='LAN')
     else:
         print 'LAN testing is off '\
           '([parallel-LAN][enable_unit_test] = 0 in ~/.pyhrf/config.cfg'

 def test_default_jde_cmd_parallel_local(self):
     t = FMRITreatment(make_outputs=False, result_dump_file=None)
     t.enable_draft_testing()
     t_fn = op.join(self.tmp_dir, 'treatment.pck')
     fout = open(t_fn, 'w')
     cPickle.dump(t, fout)
     fout.close()
     cmd = 'pyhrf_jde_estim -t %s -x local' %t_fn
     if os.system(cmd) != 0:
         raise Exception('"' + cmd + '" did not execute correctly')

def glm_analyse(fdata, contrasts, output_dir, output_prefix,
                hrf_model="Canonical", fir_delays=None,
                rescale_factor_file=None):
    glm_analyser = GLMAnalyser(hrf_model=hrf_model, contrasts=contrasts,
                               outputPrefix=output_prefix, fir_delays=fir_delays,
                               rescale_factor_file=rescale_factor_file)
    glm_analyser.set_pass_errors(False)
    glm_analyser.set_gzip_outputs(True)
    tt = FMRITreatment(fdata, glm_analyser, output_dir=output_dir)
    tt.run()

    def test_jde_estim_from_treatment_pck(self):

        t = FMRITreatment(make_outputs=False, result_dump_file=None)
        t.enable_draft_testing()
        sub_ts = t.split()
        sub_t_fn = op.join(self.tmp_dir, 'treatment.pck')
        fout = open(sub_t_fn, 'w')
        cPickle.dump(sub_ts[0], fout)
        fout.close()
        cmd = 'pyhrf_jde_estim  -t %s' %sub_t_fn
        if os.system(cmd) != 0:
            raise Exception('"' + cmd + '" did not execute correctly')

 def test_jdevemanalyser(self):
     """ Test BOLD VEM sampler on small simulation with small
     nb of iterations. Estimation accuracy is not tested.
     """
     jde_vem_analyser = JDEVEMAnalyser(beta=.8, dt=.5, hrfDuration=25.,
                                       nItMax=2, nItMin=2, fast=True,
                                       computeContrast=False, PLOT=False,
                                       constrained=True)
     tjde_vem = FMRITreatment(fmri_data=self.data_simu,
                              analyser=jde_vem_analyser,
                              output_dir=None)
     tjde_vem.run()

def jde_analyse(fdata, contrasts, output_dir):
    from pyhrf.jde.models import BOLDGibbsSampler as BG
    from pyhrf.jde.hrf import RHSampler
    from pyhrf.jde.nrl.bigaussian import NRLSampler

    sampler = BG(nb_iterations=250,
                 hrf_var=RHSampler(do_sampling=False, val_ini=np.array([0.05])),
                 response_levels=NRLSampler(contrasts=contrasts))

    analyser = JDEMCMCAnalyser(sampler=sampler)
    analyser.set_gzip_outputs(True)
    tt = FMRITreatment(fdata, analyser, output_dir=output_dir)
    tt.run(parallel='local')

 def test_jdevemanalyser(self):
     """ Test BOLD VEM sampler on small simulation with small
     nb of iterations. Estimation accuracy is not tested.
     """
     # pyhrf.verbose.set_verbosity(0)
     pyhrf.logger.setLevel(logging.WARNING)
     jde_vem_analyser = JDEVEMAnalyser(beta=.8, dt=.5, hrfDuration=25.,
                                       nItMax=2, nItMin=2, fast=True,
                                       PLOT=False,
                                       constrained=True)
     tjde_vem = FMRITreatment(fmri_data=self.data_simu,
                              analyser=jde_vem_analyser,
                              output_dir=None)
     tjde_vem.run()

    def test_default_jde_small_simulation(self):
        """ Test JDE multi-sessions sampler on small
        simulation with small nb of iterations.
        Estimation accuracy is not tested.
        """

        sampler = BMSS()

        analyser = JDEMCMCAnalyser(sampler=sampler, osfMax=4, dtMin=.4,
                                   dt=.5, driftParam=4, driftType='polynomial',
                                   outputPrefix='jde_MS_mcmc_',
                                   randomSeed=9778946)

        treatment = FMRITreatment(fmri_data=self.data_simu,
                                  analyser=analyser)

        treatment.run()

    def test_default_jde_small_simulation(self):
        """ Test ASL Physio sampler on small simulation with small nb of
        iterations. Estimation accuracy is not tested.
        """
        pyhrf.verbose.set_verbosity(0)

        sampler_params = {
            'nb_iterations' : 3,
            'smpl_hist_pace' : 1,
            'obs_hist_pace' : 1,
            'brf' : jde_asl_physio.PhysioBOLDResponseSampler(zero_constraint=False),
            'brf_var' : jde_asl_physio.PhysioBOLDResponseVarianceSampler(val_ini=\
                                                                         np.array([1e-3])),
            'prf' : jde_asl_physio.PhysioPerfResponseSampler(zero_constraint=False),
            'prf_var' : jde_asl_physio.PhysioPerfResponseVarianceSampler(val_ini=\
                                                                         np.array([1e-3])),
            'noise_var' : jde_asl_physio.NoiseVarianceSampler(),
            'drift_var' : jde_asl_physio.DriftVarianceSampler(),
            'drift' : jde_asl_physio.DriftCoeffSampler(),
            'bold_response_levels' : jde_asl_physio.BOLDResponseLevelSampler(),
            'perf_response_levels' : jde_asl_physio.PerfResponseLevelSampler(),
            'bold_mixt_params' : jde_asl_physio.BOLDMixtureSampler(),
            'perf_mixt_params' : jde_asl_physio.PerfMixtureSampler(),
            'labels' : jde_asl_physio.LabelSampler(),
            'perf_baseline' : jde_asl_physio.PerfBaselineSampler(),
            'perf_baseline_var' : jde_asl_physio.PerfBaselineVarianceSampler(),
            'check_final_value' : None,
        }


        sampler = jde_asl_physio.ASLPhysioSampler(**sampler_params)

        simu_items = phym.simulate_asl_physio_rfs(spatial_size='random_small')
        simu_fdata = FmriData.from_simulation_dict(simu_items)

        dt = simu_items['dt']
        analyser = JDEMCMCAnalyser(sampler=sampler, osfMax=4, dtMin=.4,
                                   dt=dt, driftParam=4, driftType='polynomial',
                                   outputPrefix='jde_mcmc_')

        treatment = FMRITreatment(fmri_data=simu_fdata, analyser=analyser,
                                  output_dir=None)
        treatment.run()

    def test_default_jde_small_simulation(self):
        """ Test ASL Physio sampler on small simulation with small nb of
        iterations. Estimation accuracy is not tested.
        """
        pyhrf.verbose.set_verbosity(0)

        sampler_params = {
            jde_asl_physio.ASLPhysioSampler.P_NB_ITERATIONS : 100,
            jde_asl_physio.ASLPhysioSampler.P_SMPL_HIST_PACE : 1,
            jde_asl_physio.ASLPhysioSampler.P_OBS_HIST_PACE : 1,
            'brf' : jde_asl_physio.PhysioBOLDResponseSampler(zero_constraint=False),
            'brf_var' : jde_asl_physio.PhysioBOLDResponseVarianceSampler(val_ini=\
                                                                         np.array([1e-3])),
            'prf' : jde_asl_physio.PhysioPerfResponseSampler(zero_constraint=False),
            'prf_var' : jde_asl_physio.PhysioPerfResponseVarianceSampler(val_ini=\
                                                                         np.array([1e-3])),
            'noise_var' : jde_asl_physio.NoiseVarianceSampler(),
            'drift_var' : jde_asl_physio.DriftVarianceSampler(),
            'drift_coeff' : jde_asl_physio.DriftCoeffSampler(),
            'brl' : jde_asl_physio.BOLDResponseLevelSampler(),
            'prl' : jde_asl_physio.PerfResponseLevelSampler(),
            'bold_mixt_params' : jde_asl_physio.BOLDMixtureSampler(),
            'perf_mixt_params' : jde_asl_physio.PerfMixtureSampler(),
            'label' : jde_asl_physio.LabelSampler(),
            'perf_baseline' : jde_asl_physio.PerfBaselineSampler(),
            'perf_baseline_var' : jde_asl_physio.PerfBaselineVarianceSampler(),
            'assert_final_value_close_to_true' : False,
        }


        sampler = jde_asl_physio.ASLPhysioSampler(sampler_params)

        simu_items = phym.simulate_asl_physio_rfs(spatial_size='random_small')
        simu_fdata = FmriData.from_simulation_dict(simu_items)

        dt = simu_items['dt']
        analyser = JDEMCMCAnalyser(sampler=sampler, osfMax=4, dtMin=.4,
                                   dt=dt, driftParam=4, driftType='polynomial',
                                   outputFile=None,outputPrefix='jde_mcmc_',
                                   randomSeed=None)

        treatment = FMRITreatment(fmri_data=simu_fdata, analyser=analyser)
        treatment.run()

def dummy_jde(fmri_data, dt):
    print 'run dummy_jde ...'
    jde_mcmc_sampler = \
        physio_build_jde_mcmc_sampler(3, 'basic_regularized', 
                                        do_sampling_prf=False,
                                        do_sampling_brf=False,
                                        do_sampling_prls=False,
                                        do_sampling_labels=False,
                                        do_sampling_prf_var=False,
                                        do_sampling_brf_var=False,
                                        brf_var_ini=np.array([0.1]),
                                        prf_var_ini=np.array([0.1]))

    analyser = JDEMCMCAnalyser(jde_mcmc_sampler, copy_sampler=False, dt=dt)
    analyser.set_pass_errors(False)
    tjde_mcmc = FMRITreatment(fmri_data, analyser, output_dir=None)
    outputs, fns = tjde_mcmc.run()
    print 'dummy_jde done!'
    return tjde_mcmc.analyser.sampler

 def _test_specific_parameters(self, parameter_name, fdata, simu,
                               beta=.8, dt=.5, nItMax=100, nItMin=10,
                               hrfDuration=25., estimateSigmaH=False,
                               estimateBeta=False, estimateSigmaG=False,
                               PLOT=False, constrained=True, fast=False,
                               estimateH=False, estimateG=False,
                               estimateA=False, estimateC=False,
                               estimateZ=False, estimateLA=False,
                               estimateNoise=False, estimateMP=True):
     """
     Test specific samplers.
     """
     logger.info('_test_specific_parameters %s', str(parameter_name))
     output_dir = self.tmp_dir
     # JDE analysis
     jde_vem_analyser = JDEVEMAnalyser(beta=beta, dt=dt,
                                       hrfDuration=hrfDuration,
                                       estimateSigmaH=estimateSigmaH,
                                       nItMax=nItMax, nItMin=nItMin,
                                       estimateBeta=estimateBeta,
                                       estimateSigmaG=estimateSigmaG,
                                       PLOT=PLOT,
                                       constrained=constrained, fast=fast,
                                       fmri_data=fdata,
                                       simulation=simu,
                                       estimateH=estimateH,
                                       estimateG=estimateG,
                                       estimateA=estimateA,
                                       estimateC=estimateC,
                                       estimateLabels=estimateZ,
                                       estimateLA=estimateLA,
                                       estimateMixtParam=estimateMP,
                                       estimateNoise=estimateNoise)
     tjde_vem = FMRITreatment(fmri_data=fdata, analyser=jde_vem_analyser,
                              output_dir=output_dir)
     outputs = tjde_vem.run()
     print 'out_dir:', output_dir
     return outputs

def main():
    """Run when calling the script"""

    start_time = time.time()

    if not os.path.isdir(config["output_dir"]):
        try:
            os.makedirs(config["output_dir"])
        except OSError as e:
            print("Ouput directory could not be created.\n"
                  "Error was: {}".format(e.strerror))
            sys.exit(1)

    bold_data = FmriData.from_vol_files(
        mask_file=config["parcels_file"], paradigm_csv_file=config["onsets_file"],
        bold_files=config["bold_data_file"], tr=config["tr"]
    )

    compute_contrasts, contrasts_def = load_contrasts_definitions(config["def_contrasts_file"])

    jde_vem_analyser = JDEVEMAnalyser(
        hrfDuration=config["hrf_duration"], sigmaH=config["sigma_h"], fast=True,
        computeContrast=compute_contrasts, nbClasses=2, PLOT=False,
        nItMax=config["nb_iter_max"], nItMin=config["nb_iter_min"], scale=False,
        beta=config["beta"], estimateSigmaH=True, estimateHRF=config["estimate_hrf"],
        TrueHrfFlag=False, HrfFilename='hrf.nii', estimateDrifts=True,
        hyper_prior_sigma_H=config["hrf_hyperprior"], dt=config["dt"], estimateBeta=True,
        contrasts=contrasts_def, simulation=False, estimateLabels=True,
        LabelsFilename=None, MFapprox=False, estimateMixtParam=True,
        constrained=False, InitVar=0.5, InitMean=2.0, MiniVemFlag=False, NbItMiniVem=5,
        zero_constraint=config["zero_constraint"], drifts_type=config["drifts_type"]
    )

    processing_jde_vem = FMRITreatment(
        fmri_data=bold_data, analyser=jde_vem_analyser,
        output_dir=config["output_dir"], make_outputs=True
    )

    if not config["parallel"]:
        processing_jde_vem.run()
    else:
        processing_jde_vem.run(parallel="local")


    if config["save_processing_config"]:
        # Let's canonicalize all paths
        config_save = dict(config)
        for file_nb, bold_file in enumerate(config_save["bold_data_file"]):
            config_save["bold_data_file"][file_nb] = os.path.abspath(bold_file)
        config_save["parcels_file"] = os.path.abspath(config_save["parcels_file"])
        config_save["onsets_file"] = os.path.abspath(config_save["onsets_file"])
        if config_save["def_contrasts_file"]:
            config_save["def_contrasts_file"] = os.path.abspath(config_save["def_contrasts_file"])
        config_save["output_dir"] = os.path.abspath(config_save["output_dir"])
        config_save_filename = "{}_processing.json".format(
            datetime.datetime.today()
        ).replace(" ", "_")
        config_save_path = os.path.join(config["output_dir"], config_save_filename)
        with open(config_save_path, 'w') as json_file:
            json.dump(config_save, json_file, sort_keys=True, indent=4)

    print("")
    print("Total computation took: {} seconds".format(format_duration(time.time() - start_time)))

def rfir_analyse(fdata, output_dir):
    analyser = RFIRAnalyser(RFIREstim(nb_its_max=150))
    analyser.set_gzip_outputs(True)
    tt = FMRITreatment(fdata, analyser, output_dir=output_dir)
    tt.run()

    def _test_specific_samplers(self, sampler_names, fdata,
                                nb_its=None, use_true_val=None,
                                save_history=False, check_fv=None,
                                normalize_brf=1., normalize_prf=1.,
                                normalize_mu=1., prf_prior_type='regularized',
                                brf_prior_type='regularized',
                                mu_prior_type='regularized'):
        """
        Test specific samplers.
        """
        if use_true_val is None:
            use_true_val = dict((n, False) for n in sampler_names)

        logger.info('_test_specific_samplers %s ...', str(sampler_names))

        params = deepcopy(self.sampler_params_for_single_test)

        # Loop over given samplers to enable them
        for var_name in sampler_names:
            var_class = params[var_name].__class__
            use_tval = use_true_val[var_name]

            # special case for HRF -> normalization and prior type
            if var_class == jaslh.PhysioBOLDResponseSampler:
                params[var_name] = \
                    jaslh.PhysioBOLDResponseSampler(do_sampling=True,
                                                    use_true_value=use_tval,
                                                    normalise=normalize_brf,
                                                    zero_constraint=False,
                                                    prior_type=brf_prior_type)
            elif var_class == jaslh.PhysioPerfResponseSampler:
                params[var_name] = \
                    jaslh.PhysioPerfResponseSampler(do_sampling=True,
                                                    use_true_value=use_tval,
                                                    normalise=normalize_brf,
                                                    zero_constraint=False,
                                                    prior_type=prf_prior_type)
            elif var_class == jaslh.PhysioTrueBOLDResponseSampler:
                params[var_name] = \
                    jaslh.PhysioTrueBOLDResponseSampler(do_sampling=True,
                                                        use_true_value=use_tval,
                                                        normalise=normalize_mu,
                                                        zero_constraint=False,
                                                        prior_type=mu_prior_type)

            else:
                params[var_name] = var_class(do_sampling=True,
                                             use_true_value=use_tval)

        if nb_its is not None:
            params['nb_iterations'] = nb_its

        if save_history:
            params['smpl_hist_pace'] = 1
            params['obs_hist_pace'] = 1

        if check_fv is not None:
            params['check_final_value'] = check_fv

        sampler = jaslh.ASLPhysioSampler(**params)

        output_dir = self.tmp_dir

        analyser = JDEMCMCAnalyser(sampler=sampler, osfMax=4, dtMin=.4,
                                   dt=fdata.simulation[0]['dt'], driftParam=4,
                                   driftType='polynomial',
                                   outputPrefix='jde_mcmc_',
                                   pass_error=False)

        treatment = FMRITreatment(fmri_data=fdata, analyser=analyser,
                                  output_dir=output_dir)

        outputs = treatment.run()
        print 'out_dir:', output_dir
        return outputs

    def test_default_treatment(self):

        #pyhrf.verbose.set_verbosity(4)
        t = FMRITreatment(make_outputs=False, result_dump_file=None)
        t.enable_draft_testing()
        t.run()

 def test_pickle_treatment(self):
     t = FMRITreatment(make_outputs=False, result_dump_file=None)
     t.enable_draft_testing()
     cPickle.loads(cPickle.dumps(t))

 def test_default_treatment_parallel_local(self):
     t = FMRITreatment(make_outputs=False, result_dump_file=None)
     t.enable_draft_testing()
     t.run(parallel='local')

    def _test_specific_samplers(self, sampler_names, simu,
                                nb_its=None, use_true_val=None,
                                save_history=False, check_fv=None,
                                normalize_hrf=1.,
                                hrf_prior_type='singleHRF',
                                normalize_hrf_group=1.,
                                hrf_group_prior_type='singleHRF',
                                reg_hgroup=True):
        """
        Test specific samplers.
        """
        if use_true_val is None:
            use_true_val = dict( (n,False) for n in sampler_names )

        pyhrf.verbose(1, '_test_specific_samplers %s ...' %str(sampler_names))

        params = deepcopy(self.sampler_params_for_single_test)

        # Loop over given samplers to enable them
        for var_name in sampler_names:
            var_class = params[var_name].__class__
            use_tval = use_true_val[var_name]

            # special case for HRF -> normalization and prior type
            if var_class == jms.HRF_Sampler:
                params[var_name] = jms.HRF_Sampler(do_sampling=True,
                                                   use_true_value=use_tval,
                                                   normalise=normalize_hrf,
                                                   prior_type=hrf_prior_type,
                                                   zero_contraint=False)
            elif var_class == jms.HRF_Group_Sampler:
                ptype = hrf_group_prior_type
                nhg = normalize_hrf_group
                shg = jms.HRF_Group_Sampler(do_sampling=True,
                                            use_true_value=use_tval,
                                            normalise=nhg,
                                            prior_type=ptype,
                                            zero_contraint=False,
                                            regularise=reg_hgroup)
                params[var_name] = shg
            else:
                params[var_name] = var_class(do_sampling=True,
                                             use_true_value=use_tval)


        if nb_its is not None:
            params[BMSS.P_NB_ITERATIONS] = nb_its

        if save_history:
            params[BMSS.P_SMPL_HIST_PACE] = 1
            params[BMSS.P_OBS_HIST_PACE] = 1

        if check_fv is not None:
            params[BMSS.P_CHECK_FINAL_VALUE] = check_fv

        sampler = BMSS(params)

        output_dir = self.simu_dir

        analyser = JDEMCMCAnalyser(sampler=sampler, osfMax=4, dtMin=.4,
                                   dt=.6, driftParam=4, driftType='polynomial',
                                   outputFile=None, outputPrefix='jde_mcmc_',
                                   randomSeed=5421087, pass_error=False)

        treatment = FMRITreatment(fmri_data=simu, analyser=analyser,
                                  output_dir=output_dir)


        outputs = treatment.run()
        #print 'out_dir:', output_dir
        return outputs