def getOutputs(self):
        outputs = {}
        logger.info('get output of sampled variables ...')
        for v in self.variables:
            logger.debug('get outputs of %s', v.name)
            outputs.update(v.getOutputs())

        logger.info('get output of global observables ...')
        outputs.update(self.getGlobalOutputs())

        #tp = time.time()
        d = {'parcel_size': np.array([self.dataInput.nbVoxels])}
        outputs['analysis_duration'] = xndarray(\
                                            np.array([self.analysis_duration]),
                                            axes_names=['parcel_size'],
                                            axes_domains=d)
        try:
            outputs['conv_error'] = xndarray(np.array(self.converror))
            outputs['loglikelihood'] = xndarray(np.array([self.loglikelihood]))
            outputs['bic'] = xndarray(np.array([self.bic]),
                                                axes_names = ['parcel_size'],
                                                axes_domains = d)
        except AttributeError:
            pass

        # for on, o in outputs.iteritems():
        #     #print 'on:', o
        #     yield (on,o)
        return outputs

    def test_stack(self):

        d1 = np.arange(4 * 5).reshape(4, 5)
        c1 = xndarray(d1, ['x', 'y'],
                      {'x': np.arange(4) * 3,
                       'y': np.arange(5) * 3,
                       })
        if debug:
            print 'c1:'
            print c1.descrip()

        d2 = np.arange(4 * 5).reshape(4, 5) * 2
        c2 = xndarray(d2, ['x', 'y'],
                      {'x': np.arange(4) * 3,
                       'y': np.arange(5) * 3,
                       })
        if debug:
            print 'c2:'
            print c2.descrip()

        c_stacked = stack_cuboids([c1, c2], 'stack_axis', ['c1', 'c2'])

        if debug:
            print 'c_stacked:'
            print c_stacked.descrip()

        a_stacked = np.array([d1, d2])
        assert (c_stacked.data == a_stacked).all()

    def _get_outputs(self, output_type='ndarray'):
        """
        output_type : 'ndarray' or 'cuboid'
        """

        outputs = {}

        # Default outputs: post. mean, post. var, trajectories.
        on = self.name + '_obs_mean'
        if output_type == 'cuboid':
            outputs[on] = xndarray(self.obs_mean, axes_names=self.axes_names,
                                   axes_domains=self.axes_domains)
        else:
            outputs[on] = self.obs_mean

        on = self.name + '_obs_var'
        if output_type == 'cuboid':
            outputs[on] = xndarray(self.obs_var, axes_names=self.axes_names,
                                   axes_domains=self.axes_domains)
        else:
            outputs[on] = self.obs_var

        # Custom outputs:
        self.set_outputs(outputs, output_type)

        return outputs

 def test_xmapping_inconsistent_mapping_value(self):
     a = xndarray(np.arange(2 * 4).reshape(2, 4).T, ['time', 'parcel'],
                  {'time': np.arange(4) * .5, 'parcel': [2, 6]})
     parcel_map = xndarray(np.array([[2, 2, 2, 6], [6, 6, 6, 0], [6, 6, 0, 0]]),
                           ['axial', 'coronal'], value_label='ROI')
     self.assertRaisesRegexp(ArrayMappingError,
                             r'Value label "ROI" of xmapping not found '
                             'in array axes \(time, parcel\)',
                             a.map_onto, parcel_map)

 def test_xmapping_inconsistent_domain(self):
     a = xndarray(np.arange(2 * 4).reshape(2, 4).T, ['time', 'parcel'],
                  {'time': np.arange(4) * .5, 'parcel': [2, 7]})
     parcel_map = xndarray(np.array([[2, 2, 2, 6], [6, 6, 6, 0], [6, 6, 0, 0]]),
                           ['axial', 'coronal'], value_label='parcel')
     self.assertRaisesRegexp(ArrayMappingError,
                             'Domain of axis "parcel" to be mapped is '
                             'not a subset of values in the mapping array.',
                             a.map_onto, parcel_map)

    def test_equality(self):
        c1 = xndarray(self.arr3d, self.arr3dNames, self.arr3dDom,
                      self.arr3dLabel)

        c2 = xndarray(self.arr3d, self.arr3dNames, self.arr3dDom,
                      self.arr3dLabel)

        assert c1 == c2

        c3 = xndarray(self.arr3d, self.arr3dNames)

        assert c1 != c3

def compute_T_Pvalue(betas, stds_beta, mask_file, null_hyp=True):
    '''
    Compute Tvalues statistic and Pvalue based upon estimates
    and their standard deviation
    beta and std_beta for all voxels
    beta: shape (nb_vox, 1)
    std: shape (1)
    Assume null hypothesis if null_hyp is True
    '''
    from pyhrf.ndarray import xndarray

    import sys
    sys.path.append("/home/i2bm/BrainVisa/source/pyhrf/pyhrf-free/trunk/script/WIP/Scripts_IRMf_Adultes_Solv/Scripts_divers_utiles/Scripts_utiles/")
    from Functions_fit import Permutation_test, stat_mean, stat_Tvalue, stat_Wilcoxon

    mask = xndarray.load(mask_file).data #to save P and Tval on a map

    BvalC = xndarray(betas, axes_names=['sagittal', 'coronal', 'axial'])
    Betasval = BvalC.flatten(mask, axes=['sagittal', 'coronal', 'axial'], new_axis='position').data

    Stdsval = stds_beta

    Tval = xndarray(Betasval/Stdsval, axes_names=['position']).data

    nb_vox = Betasval.shape[0]
    nb_reg = betas.shape[1]
    dof = nb_vox - nb_reg #degrees of freedom for STudent distribution
    assert dof>0

    Probas=np.zeros(Betasval.shape)
    for i in xrange(nb_vox):
        if null_hyp:
            #STudent distribution
            from scipy.stats import t
            fmix = lambda x: t.pdf(x, dof)
        else:
            fmix = lambda t:  1/np.sqrt(2*np.pi*Stdsval[i]**2)*np.exp(- (t - Betasval[i])**2 / (2*Stdsval[i]**2) )
        Probas[i] = quad(fmix, Tval[i], float('inf'))[0]

    Tvalues_ = xndarray(Tval, axes_names=['position'])
    Pvalues_ = xndarray(Probas, axes_names=['position'])
    Tvalues = Tvalues_.expand(mask, 'position', ['sagittal','coronal','axial'])
    Pvalues = Pvalues_.expand(mask, 'position', ['sagittal','coronal','axial'])

    #Computation of Pvalue using permutations
    #not possible to do this actually...it was used for group level stats
    #Pvalue_t = np.zeros(Betasval.shape)
    #for i in xrange(nb_vox):
        #Pvalue_t[i] = Permutation_test(Betasval[i], n_permutations=10000, \
                    #stat = stat_Tvalue, two_tailed=False, plot_histo=False)

    return Tvalues.data, Pvalues.data

    def test_split(self):

        # pyhrf.verbose.set_verbosity(0)
        pyhrf.logger.setLevel(logging.WARNING)
        sh = (2, 4, 4, 4)
        c = xndarray(np.arange(np.prod(sh)).reshape(sh), ['condition'] + MRI3Daxes,
                     {'condition': ['audio', 'video']})
        if debug:
            print 'Original cub:'
            print c.descrip()

        fn = op.join(self.tmp_dir, 'cub.nii')
        if debug:
            print 'Save and load original cub'
        c.save(fn)
        c = xndarray.load(fn)

        fn = op.join(self.tmp_dir, 'cub2.nii')
        if debug:
            print 'Save and load new cub with meta data from original cuboid'
        sh = (4, 4, 4)
        c2 = xndarray(np.arange(np.prod(sh)).reshape(sh), MRI3Daxes,
                      meta_data=c.meta_data)
        c2.save(fn)
        c2 = xndarray.load(fn)

        fns = []
        sub_cuboids = []
        if debug:
            print 'Split and save sub cuboids'
        for dvalue, sub_c in c.split('condition').iteritems():
            fn = op.join(self.tmp_dir, add_suffix('sub_c.nii',
                                                  '_%s' % str(dvalue)))
            if debug and dvalue == 'audio':

                print 'fn_out:', fn
                print 'sub_c:'
                print sub_c.descrip()
                sub_cuboids.append(sub_c)
                sub_c.save(fn)
                fns.append(fn)
        if debug:
            print ''
            print 'Load sub c again ...'
        for fn, sub_c in zip(fns, sub_cuboids):
            if debug:
                print 'fn:', fn
            c_loaded = xndarray.load(fn)
            if debug:
                print 'c_loaded:'
                print c_loaded.descrip()
            self.assertEqual(c_loaded, sub_c)

    def test_xmapping(self):
        a = xndarray(np.arange(2 * 4).reshape(2, 4).T, ['time', 'parcel'],
                     {'time': np.arange(4) * .5, 'parcel': [2, 6]})
        parcel_map = xndarray(np.array([[2, 2, 2, 6], [6, 6, 6, 0], [6, 6, 0, 0]]),
                              ['axial', 'coronal'], value_label='parcel')
        a_mapped = a.map_onto(parcel_map)

        self.assertEqual(a_mapped.data.shape,
                         parcel_map.data.shape + a.data.shape[:1])
        self.assertEqual(a_mapped.axes_names, ['axial', 'coronal', 'time'])
        npt.assert_array_equal(a_mapped.get_domain('time'),
                               a.get_domain('time'))
        npt.assert_array_equal(a_mapped.data[0, 0], a.data[:, 0])
        npt.assert_array_equal(a_mapped.data[1, 0], a.data[:, 1])
        npt.assert_array_equal(a_mapped.data[-1, -1], 0.)

def view(data, axesNames=None, axesDomains=None, errors=None,
         valueLabel='value', mask=None, maskAxes=None, maskName='mask'):

    """
    Interactively browse and display a n-dimensional numpy array. Axes can be
    labeled with a name in 'axesNames' and associated to real values in
    'axesDomains'. Errors data with the same data type and shape as 'data' can
    be given with 'errors' argument.
    Example :
    from numpy import *
    import ndview
    a = array([ [4,5,6],[8,10,12] ])
    names = ['time','position']
    domains = {'time':[0.1, 0.2]}
    ndview.view(a, axesNames=names, axesDomains=domains)
    """
    if data.dtype == np.int16:
        data = data.astype(np.int32)
    if not isinstance(data, xndarray):
        c = xndarray(data, errors=errors, axesDomains=axesDomains,
                   axesNames=axesNames,
                   value_label=valueLabel)
    else:
        c = data
    viewxndarray(c, mask=mask, maskName=maskName, maskAxes=maskAxes)

    def test_save_as_nii(self):

        c = xndarray(self.arr3d, ['x', 'y', 'z'],
                     {'x': np.arange(self.sh3d[0]) * 3,
                      'y': np.arange(self.sh3d[1]) * 3,
                      'z': np.arange(self.sh3d[2]) * 3,
                      },
                     self.arr3dLabel)

        if debug:
            print 'Original cuboid:'
            print c.descrip()

        fn = op.join(self.tmp_dir, 'cuboid.nii')
        if debug:
            print 'fn:', fn

        c.save(fn)
        assert op.exists(fn)

        c_loaded = xndarray.load(fn)

        if debug:
            print 'Loaded cuboid:'
            print c_loaded.descrip()

        assert c == c_loaded

    def test_plot_cuboid_as_curve(self):
        from pyhrf.ndarray import xndarray

        sh = (10, 10, 5, 3)
        data = np.zeros(sh)
        data[:, :, :, 0] = 1.0
        data[:, :, :, 1] = 2.0
        data[:, :, :, 2] = 3.0
        c1 = xndarray(
            data,
            axes_names=["sagittal", "coronal", "axial", "condition"],
            axes_domains={"condition": ["audio1", "audio2", "video"]},
        )

        f = plt.figure()
        ax = f.add_subplot(111)

        ori = ["condition", "sagittal"]
        plot_cub_as_curve(
            c1.sub_cuboid(axial=0, coronal=0).reorient(ori),
            colors={"audio1": "red", "audio2": "orange", "video": "blue"},
            axes=ax,
        )
        if 0:
            plt.show()

    def test_plot_cuboid2d_as_image(self):
        from pyhrf.ndarray import xndarray
        import matplotlib

        sh = (10, 3)
        c1 = xndarray(
            np.arange(np.prod(sh)).reshape(sh),
            axes_names=["sagittal", "condition"],
            axes_domains={"condition": ["audio1", "audio2", "video"]},
        )

        f = plt.figure()
        ax = f.add_subplot(111)

        ori = ["condition", "sagittal"]
        cm = matplotlib.cm.get_cmap("winter")
        norm = matplotlib.colors.Normalize(vmin=5, vmax=20)
        plot_cub_as_image(
            c1.reorient(ori),
            cmap=cm,
            norm=norm,
            axes=ax,
            show_axes=True,
            show_axis_labels=True,
            show_colorbar=True,
            show_tick_labels=False,
        )
        if 0:
            plt.show()

    def test_to_latex_3d_hide_name_style(self):

        sh = (2, 3, 4)
        c = xndarray(np.arange(np.prod(sh)).reshape(sh), ['mod', 'time', 'subj'],
                     {'mod': ['m1', 'm2'],
                      'time': np.arange(3) * .5,
                      'subj': ['s1', 's2', 's3', 's4']})

        if debug:
            print 'c:'
            print c.descrip()

        result = c.to_latex(col_axes=['mod', 'subj'], row_axes=['time'],
                            header_styles={'mod': 'join', 'time': 'hide_name'},
                            hval_fmt={'time': '%1.1f'}, val_fmt='%d')

        l1 = range(4) + range(12, 16)
        l2 = range(4, 8) + range(16, 20)
        l3 = range(8, 12) + range(20, 24)
        good_result = \
            '\\begin{tabular}{c | ' + ' '.join(['c'] * 8) + '}\n' \
            '&\multicolumn{4}{c}{mod=m1} & \multicolumn{4}{c}{mod=m2}\\\\\n' \
            '&\multicolumn{8}{c}{subj}\\\\\n' \
            '&s1 & s2 & s3 & s4 & s1 & s2 & s3 & s4\\\\\n' + \
            '0.0 & ' + \
            ' & '.join([str(i) for i in l1]) + '\\\\\n' + \
            '0.5 & ' + ' & '.join([str(i) for i in l2]) + '\\\\\n' + \
            '1.0 & ' + ' & '.join([str(i) for i in l3]) + '\\\\\n' + \
            '\\end{tabular}'

        self.assertEqual(result, good_result,
                         'to latex returned:\n' + result + '\n' + 'expected:\n' + good_result)

    def test_to_latex_3d_inner_axes(self):

        sh = (2, 3, 4)
        c = xndarray(np.arange(np.prod(sh)).reshape(sh), ['mod', 'time', 'subj'],
                     {'mod': ['m1', 'm2'],
                      'time': np.arange(3) * .5,
                      'subj': ['s1', 's2', 's3', 's4']})

        if debug:
            print 'c:'
            print c.descrip()

        result = c.to_latex(col_axes=['subj'], row_axes=['time'],
                            inner_axes=['mod'],
                            header_styles={'mod': 'join', 'time': 'join'},
                            val_fmt="%d",
                            hval_fmt={'time': "%1.1f"})

        l1 = ['%d | %d' % (a, b) for a, b in zip(range(4), range(12, 16))]
        l2 = ['%d | %d' % (a, b) for a, b in zip(range(4, 8), range(16, 20))]
        l3 = ['%d | %d' % (a, b) for a, b in zip(range(8, 12), range(20, 24))]
        good_result = \
            '\\begin{tabular}{c | ' + ' '.join(['c'] * 4) + '}\n' \
            '&\multicolumn{4}{c}{subj}\\\\\n' \
            '&s1 & s2 & s3 & s4\\\\\n' + \
            'time=0.0 & ' + \
            ' & '.join([str(i) for i in l1]) + '\\\\\n' + \
            'time=0.5 & ' + ' & '.join([str(i) for i in l2]) + '\\\\\n' + \
            'time=1.0 & ' + ' & '.join([str(i) for i in l3]) + '\\\\\n' + \
            '\\end{tabular}'

        self.assertEqual(result, good_result,
                         'to latex returned:\n' + result + '\n' + 'expected:\n' + good_result)

    def test_save_as_gii(self):

        c = xndarray(np.arange(4 * 2).reshape(4, 2), ['x', 'time'],
                     {'time': np.array([2 / 3., 2.3])})

        if debug:
            print 'Original cuboid:'
            print c.descrip()
            print c.data

        fn = op.join(self.tmp_dir, 'cuboid.gii')
        if debug:
            print 'fn:', fn

        c.save(fn)
        assert op.exists(fn)

        c_loaded = xndarray.load(fn)

        if debug:
            print 'Loaded cuboid:'
            print c_loaded.descrip()
            print c_loaded.data

        assert c == c_loaded

    def test_save_as_nii(self):

        # pyhrf.verbose.set_verbosity(0)
        pyhrf.logger.setLevel(logging.WARNING)
        c = xndarray(self.arr3d, ['x', 'y', 'z'],
                     {'x': np.arange(self.sh3d[0]) * 3,
                      'y': np.arange(self.sh3d[1]) * 3,
                      'z': np.arange(self.sh3d[2]) * 3,
                      },
                     self.arr3dLabel)

        if debug:
            print 'Original cuboid:'
            print c.descrip()

        fn = op.join(self.tmp_dir, 'cuboid.nii')
        if debug:
            print 'fn:', fn

        c.save(fn)
        assert op.exists(fn)

        c_loaded = xndarray.load(fn)

        if debug:
            print 'Loaded cuboid:'
            print c_loaded.descrip()

        assert c == c_loaded

    def test_operations(self):

        data1 = np.arange(4 * 5, dtype=float).reshape(4, 5) + 1
        data1_copy = data1.copy()
        data2 = np.arange(4 * 5).reshape(4, 5) * 4. + 1
        data2_copy = data2.copy()

        c = xndarray(data1, ['x', 'time'], {'time': np.arange(5) + 2 / 3.})
        c_bak = c

        c2 = xndarray(data2, ['x', 'time'], {'time': np.arange(5) + 2 / 3.})

        c_add = c + c2
        assert (c_add.data == data1_copy + data2_copy).all()

        c += c2
        assert c.data is data1
        assert (c.data == data1_copy + data2_copy).all()
        assert c is c_bak

        c_sub = c - c2
        assert (c_sub.data == data1 - data2).all()

        c -= c2
        assert (c_sub.data == data1_copy).all()
        assert c is c_bak

        c_mul = c * c2
        assert (c_mul.data == data1_copy * data2_copy).all()

        c_div = c / c2
        assert (c_div.data.flat[1:] == (
            data1_copy / data2_copy).flat[1:]).all()

        r = 2 / c
        npt.assert_array_equal(r.data, 2 / c.data)

        r = 4 * c
        npt.assert_array_equal(r.data, 4 * c.data)

        r = 4 + c
        npt.assert_array_equal(r.data, 4 + c.data)

        r = 5 - c
        npt.assert_array_equal(r.data, 5 - c.data)

    def setUp(self):
        self.cub0 = xndarray(_np.random.rand(10,10))
        self.cub3DVol = xndarray(_np.random.rand(10,10,10),
                               axes_names=MRI3Daxes)
        d4D = _np.zeros((2,2,2,3))
        for c in xrange(3):
            d4D[:,:,:,c] = _np.ones((2,2,2))*(c-2)

        self.cub4DVol = xndarray(d4D, axes_names=['condition']+MRI3Daxes)

        self.cub4DTimeVol = xndarray(_np.random.rand(100,10,10,10),
                               axes_names=['time']+MRI3Daxes)
        self.cubNDVol = xndarray(_np.random.rand(10,2,2,2,3),
                               axes_names=['time']+MRI3Daxes+['condition'],
                               axes_domains={'condition':['audio','video','na']})

        self.tmp_dir = tempfile.mkdtemp(prefix='pyhrf_tests',
                                        dir=pyhrf.cfg['global']['tmp_path'])

    def test_unstack_empty_inner_axes(self):
        size = 6
        ad = ['out_dv%d'%i for i in range(size)]
        c = xndarray(np.arange(size), axes_names=['a1'], axes_domains={'a1':ad})
        uc = c.unstack(['a1'], [])

        self.assertEqual(uc.data.shape, (size,))
        self.assertEqual(len(uc.axes_domains), 1)
        npt.assert_array_equal(uc.axes_domains['a1'], ad)