def func(op1, op2):
        op = op1(op2)
        attr = {}
        attr.update(op2.attrout)
        attr.update(op1.attrout)
        assert_equal(op.attrout, attr)
        assert_is(op.classout, op1.classout)
        if op1.flags.linear:
            assert_is_type(op, ZeroOperator)
            assert_same(op.todense(shapein=3, shapeout=4), np.zeros((4, 3)))
            return
        if op1.flags.shape_output == 'unconstrained' or \
           op1.flags.shape_input != 'explicit' and \
           op2.flags.shape_output != 'explicit':
            assert_is_type(op, CompositionOperator)
        else:
            assert_is_type(op, ConstantOperator)

        if op1.flags.shape_input == 'unconstrained' and \
           op2.flags.shape_output == 'unconstrained':
            return
        with rule_manager(none=True):
            op_ref = op1(op2)
        assert_same(op.todense(shapein=3, shapeout=4),
                    op_ref.todense(shapein=3, shapeout=4))

 def func(d):
     x_ = np.array(x if d.kind == 'c' else x.real, dtype=d)
     actual = abs2(x_)
     expected = np.abs(x_**2)
     assert_same(actual, expected)
     abs2(x_, actual)
     assert_same(actual, expected)

def test_endian():
    bendian = _read_filters(path + 'invntt_be', bigendian=True)
    lendian = _read_filters(path + 'invntt_le')
    assert len(bendian) == 6
    assert bendian[0]['data'].size == 101
    assert_same(bendian[0]['data'][0], 5597147.4155586753)
    assert_eq(lendian, bendian)

def test_sky2():
    assert_allclose(outputs['I'], outputs['IQU'][..., 0], atol=2e-2)
    assert_allclose(outputs['QU'], outputs['IQU'][..., 1:], atol=2e-2)

    for k in ('QU', 'IQU'):
        assert_equal(cbiks[k], cbiks['I'])
        assert_same(pT1s[k], pT1s['I'].astype(np.float32), rtol=15)

def test_add_subtract_grid_operator():
    acq = QubicAcquisition(150, sampling, detector_ngrids=2)
    tod = map2tod(acq, input_map, convolution=False)
    add = acq.get_add_grids_operator()
    sub = acq.get_subtract_grids_operator()
    assert_same(add(tod), tod[:992] + tod[992:])
    assert_same(sub(tod), tod[:992] - tod[992:])

def test_equ2gal():
    equ2gal = CartesianEquatorial2GalacticOperator()
    gal2equ = CartesianGalactic2EquatorialOperator()

    assert equ2gal.I is equ2gal.T
    assert gal2equ.I is gal2equ.T
    assert_same(equ2gal.todense(shapein=3), equ2gal.data)
    assert_same(gal2equ.todense(shapein=3), equ2gal.data.T)

    shapes = (3, (2, 3), (4, 2, 3))

    def func(op, shape):
        vec = np.arange(product(shape)).reshape(shape)
        vec_ = vec.reshape((-1, 3))
        expected = np.empty(shape)
        expected_ = expected.reshape((-1, 3))
        for i in range(expected.size // 3):
            expected_[i] = op(vec_[i])
        actual = op(vec)
        assert_same(actual, expected)

    for op in (equ2gal, gal2equ):
        for shape in shapes:
            yield func, op, shape

    assert_is_instance(equ2gal(gal2equ), IdentityOperator)
    assert_is_instance(gal2equ(equ2gal), IdentityOperator)

 def func(dtype, shape):
     sqrt3 = np.sqrt(np.array(3, dtype))
     c = [sqrt3 / 2, np.array(1, dtype) / 2]
     coords = np.empty(shape + (2,), dtype)
     coords[...] = c
     assert_same(rotate(coords, 30), e, broadcasting=True)
     assert_same(rotate(coords, 30, out=coords), e, broadcasting=True)

 def func(n, m):
     slices = split(n, m)
     assert_eq(len(slices), m)
     x = np.zeros(n, int)
     for s in slices:
         x[s] += 1
     assert_same(x, 1, broadcasting=True)
     assert_eq([split(n, m, i) for i in range(m)], slices)

 def func(spacing, angle):
     origin = shape_grid * np.array(spacing) * np.sqrt(2) / 2
     grid = LayoutGridSquares(shape_grid, spacing, origin=origin, angle=angle, filling_factor=filling_factor)
     shape_plane = np.ceil(shape_grid * np.array(spacing) * np.sqrt(2))
     scene = SceneGrid(shape_plane)
     proj = scene.get_integration_operator(grid.vertex)
     assert not proj.outside
     assert_same(np.sum(proj.matrix.data.value), len(grid) * spacing ** 2 * filling_factor, rtol=100)

 def func2(itype, ftype, shapein):
     mat = get_mat(itype, ftype)
     op = SparseOperator(mat)
     todense = op.todense(shapein=(6,) + shapein)
     assert_same(todense.T, op.T.todense(shapeout=(6,) + shapein))
     op2 = SparseOperator(mat, shapeout=(2, 2) + shapein,
                          shapein=(3, 2) + shapein)
     assert_same(op2.todense(), todense)

def test_scatter():
    np.random.seed(0)
    n = 4
    x = np.random.random(n)
    layout = PackedTable(n, x=x)
    s = split(n, size, rank)
    scattered = layout.scatter()
    assert_same(scattered.x, x[s])
    assert_same(scattered.all.x, x)

 def func(lmax, n):
     class XpolDummy(Xpol):
         def __init__(self):
             self.lmax = lmax
             self.wl = np.ones(max(n+1, 1))
     xpol = XpolDummy()
     mll = xpol._get_Mll(binning=False)
     expected = FitsArray('test/data/xpol_mll_{}_{}.fits'.format(lmax, n))
     assert_same(mll, expected, atol=200)

 def func(s, v):
     if s is setattr_unpacked:
         assert_raises(TypeError, s, layout, 'key', v)
         return
     s(layout, 'key', v)
     assert isalias(layout.key, val)
     assert isalias(layout.all.key, val)
     assert_same(layout.key, val)
     assert_same(layout.all.key, val.reshape(shape))

 def func(op, shape):
     vec = np.arange(product(shape)).reshape(shape)
     vec_ = vec.reshape((-1, 3))
     expected = np.empty(shape)
     expected_ = expected.reshape((-1, 3))
     for i in range(expected.size // 3):
         expected_[i] = op(vec_[i])
     actual = op(vec)
     assert_same(actual, expected)

def test_special_attribute_func3():
    shape = (6, 6)
    ordering = np.arange(product(shape))[::-1].reshape(shape)
    val = np.arange(product(shape)).reshape(shape)
    func = lambda s, x: x * s.val * s.myscalar1 * s._myscalar2
    layout = PackedTable(shape, ordering=ordering, val=val, key=func)
    layout.myscalar1 = 2
    layout._myscalar2 = 3
    assert_same(layout.key(2), val.ravel()[::-1] * 12)
    assert_same(layout.all.key(2), val * 12)

 def func(format):
     cls = getattr(sp, format + '_matrix')
     shapein = (2, 2)
     shapeout = (1, 4, 1)
     so = SparseOperator(cls(A), shapein=shapein, shapeout=shapeout)
     for vec in vecs:
         assert_same(so(np.reshape(vec, shapein)),
                     np.dot(A, vec).reshape(shapeout))
         assert_same(so.T(np.reshape(vec, shapeout)),
                     np.dot(A.T, vec).reshape(shapein))

def test_noiseless():
    sampling = create_random_pointings([0, 90], 100, 10)
    acq_qubic = QubicAcquisition(150, sampling)
    acq_planck = PlanckAcquisition(150, acq_qubic.scene, true_sky=SKY)
    acq_fusion = QubicPlanckAcquisition(acq_qubic, acq_planck)
    np.random.seed(0)
    y1 = acq_fusion.get_observation()
    np.random.seed(0)
    y2 = acq_fusion.get_observation(noiseless=True) + acq_fusion.get_noise()
    assert_same(y1, y2)

def test_gaussian():
    shape = (100, 90)
    sigma = 2.1, 2.4
    center = 14.3, 52.2
    dtype = float
    y, x = np.ogrid[0:shape[0], 0:shape[1]]
    ref = np.exp(-(x - center[0])**2 / (2 * sigma[0]**2) +
                 -(y - center[1])**2 / (2 * sigma[1]**2)) / \
        (2 * np.pi * sigma[0] * sigma[1])
    actual = gaussian(shape, sigma=sigma, center=center, dtype=dtype)
    assert_same(actual, ref, atol=2)

 def func(algo, reuse_initial_state):
     g = algo(reuse_initial_state=reuse_initial_state)
     pi = g.run()
     assert g.niterations == 10
     assert_same(pi, np.pi)
     if g.reuse_initial_state:
         assert_raises(RuntimeError, g.restart)
         return
     g.restart()
     assert g.niterations == 10
     assert_same(pi, np.pi)

 def func(shape, degrees):
     angle = np.arange(product(shape)).reshape(shape)
     if degrees:
         angle_ = np.radians(angle)
     else:
         angle_ = angle
     angle_ = angle_.reshape(angle.size)
     r = Rotation2dOperator(angle, degrees=degrees)
     actual = r([1, 0]).reshape((angle.size, 2))
     expected = np.array([np.cos(angle_), np.sin(angle_)]).T
     assert_same(actual, expected)