def test_cache_key(self):
        def fn(x):
            return x ** 2

        f = Gradient(fn)
        self.assertTrue(np.allclose(f(3), 0.0))
        self.assertTrue(np.allclose(f(3.0), 6.0))

def order_segments(segments):
    '''Piece the segments together in order, return a list of vertices
    '''
    segments = [np.array(seg).tolist() for seg in segments]
    if not segments:
        return
    verts = [segments[0][0], segments[0][1]]
    original = range(1, len(segments))
    while original:
        match = False
        for ind, segment in enumerate(segments):
            if not ind in original:
                continue
            pt1, pt2 = segment
            if np.allclose(pt1, verts[-1]):
                verts.append(pt2)
                original.remove(ind)
                match = True
            elif np.allclose(pt2, verts[-1]):
                verts.append(pt1)
                original.remove(ind)
                match = True
        if not match:
            verts.append(verts[0])
            return np.array(verts)
    return np.array(verts)

def _validate_covars(covars, cvtype, nmix, n_dim):
    from scipy import linalg
    if cvtype == 'spherical':
        if len(covars) != nmix:
            raise ValueError("'spherical' covars must have length nmix")
        elif np.any(covars <= 0):
            raise ValueError("'spherical' covars must be non-negative")
    elif cvtype == 'tied':
        if covars.shape != (n_dim, n_dim):
            raise ValueError("'tied' covars must have shape (n_dim, n_dim)")
        elif (not np.allclose(covars, covars.T)
              or np.any(linalg.eigvalsh(covars) <= 0)):
            raise ValueError("'tied' covars must be symmetric, "
                             "positive-definite")
    elif cvtype == 'diag':
        if covars.shape != (nmix, n_dim):
            raise ValueError("'diag' covars must have shape (nmix, n_dim)")
        elif np.any(covars <= 0):
            raise ValueError("'diag' covars must be non-negative")
    elif cvtype == 'full':
        if covars.shape != (nmix, n_dim, n_dim):
            raise ValueError("'full' covars must have shape "
                             "(nmix, n_dim, n_dim)")
        for n, cv in enumerate(covars):
            if (not np.allclose(cv, cv.T)
                or np.any(linalg.eigvalsh(cv) <= 0)):
                raise ValueError("component %d of 'full' covars must be "
                                 "symmetric, positive-definite" % n)

def test_massivenu_density():
    # Testing neutrino density calculation

    # Simple test cosmology, where we compare rho_nu and rho_gamma
    # against the exact formula (eq 24/25 of Komatsu et al. 2011)
    # computed using Mathematica.  The approximation we use for f(y)
    # is only good to ~ 0.5% (with some redshift dependence), so that's
    # what we test to.
    ztest = np.array([0.0, 1.0, 2.0, 10.0, 1000.0])
    nuprefac = 7.0 / 8.0 * (4.0 / 11.0) ** (4.0 / 3.0)
    #  First try 3 massive neutrinos, all 100 eV -- note this is a universe
    #  seriously dominated by neutrinos!
    tcos = core.FlatLambdaCDM(75.0, 0.25, Tcmb0=3.0, Neff=3,
                              m_nu=u.Quantity(100.0, u.eV))
    assert tcos.has_massive_nu
    assert tcos.Neff == 3
    nurel_exp = nuprefac * tcos.Neff * np.array([171969, 85984.5, 57323,
                                                 15633.5, 171.801])
    assert np.allclose(tcos.nu_relative_density(ztest), nurel_exp, rtol=5e-3)
    assert np.allclose(tcos.efunc([0.0, 1.0]), [1.0, 7.46144727668], rtol=5e-3)

    # Next, slightly less massive
    tcos = core.FlatLambdaCDM(75.0, 0.25, Tcmb0=3.0, Neff=3,
                              m_nu=u.Quantity(0.25, u.eV))
    nurel_exp = nuprefac * tcos.Neff * np.array([429.924, 214.964, 143.312,
                                                 39.1005, 1.11086])
    assert np.allclose(tcos.nu_relative_density(ztest), nurel_exp,
                       rtol=5e-3)

    # For this one also test Onu directly
    onu_exp = np.array([0.01890217, 0.05244681, 0.0638236,
                        0.06999286, 0.1344951])
    assert np.allclose(tcos.Onu(ztest), onu_exp, rtol=5e-3)

    # And fairly light
    tcos = core.FlatLambdaCDM(80.0, 0.30, Tcmb0=3.0, Neff=3,
                              m_nu=u.Quantity(0.01, u.eV))

    nurel_exp = nuprefac * tcos.Neff * np.array([17.2347, 8.67345, 5.84348,
                                                 1.90671, 1.00021])
    assert np.allclose(tcos.nu_relative_density(ztest), nurel_exp,
                       rtol=5e-3)
    onu_exp = np.array([0.00066599, 0.00172677, 0.0020732,
                        0.00268404, 0.0978313])
    assert np.allclose(tcos.Onu(ztest), onu_exp, rtol=5e-3)
    assert np.allclose(tcos.efunc([1.0, 2.0]), [1.76225893, 2.97022048],
                       rtol=1e-4)
    assert np.allclose(tcos.inv_efunc([1.0, 2.0]), [0.5674535, 0.33667534],
                       rtol=1e-4)

    # Now a mixture of neutrino masses, with non-integer Neff
    tcos = core.FlatLambdaCDM(80.0, 0.30, Tcmb0=3.0, Neff=3.04,
                              m_nu=u.Quantity([0.0, 0.01, 0.25], u.eV))
    nurel_exp = nuprefac * tcos.Neff * np.array([149.386233, 74.87915, 50.0518,
                                                 14.002403, 1.03702333])
    assert np.allclose(tcos.nu_relative_density(ztest), nurel_exp,
                       rtol=5e-3)
    onu_exp = np.array([0.00584959, 0.01493142, 0.01772291,
                        0.01963451, 0.10227728])
    assert np.allclose(tcos.Onu(ztest), onu_exp, rtol=5e-3)

    def testFull(self, num_best=None, shardsize=100):
        if self.cls == similarities.Similarity:
            index = self.cls(None, corpus, num_features=len(dictionary), shardsize=shardsize)
        else:
            index = self.cls(corpus, num_features=len(dictionary))
        if isinstance(index, similarities.MatrixSimilarity):
            expected = numpy.array([
                [0.57735026, 0.57735026, 0.57735026, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                [0.0, 0.40824831, 0.0, 0.40824831, 0.40824831, 0.40824831, 0.40824831, 0.40824831, 0.0, 0.0, 0.0, 0.0],
                [0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0],
                [0.0, 0.0, 0.40824831, 0.0, 0.0, 0.0, 0.81649661, 0.0, 0.40824831, 0.0, 0.0, 0.0],
                [0.0, 0.0, 0.0, 0.57735026, 0.57735026, 0.0, 0.0, 0.57735026, 0.0, 0.0, 0.0, 0.0],
                [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1., 0.0, 0.0],
                [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.70710677, 0.70710677, 0.0],
                [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.57735026, 0.57735026, 0.57735026],
                [0.0, 0.0, 0.0, 0.0, 0.0, 0.57735026, 0.0, 0.0, 0.0, 0.0, 0.57735026, 0.57735026],
                ], dtype=numpy.float32)
            # HACK: dictionary can be in different order, so compare in sorted order
            self.assertTrue(numpy.allclose(sorted(expected.flat), sorted(index.index.flat)))
        index.num_best = num_best
        query = corpus[0]
        sims = index[query]
        expected = [(0, 0.99999994), (2, 0.28867513), (3, 0.23570226), (1, 0.23570226)][ : num_best]

        # convert sims to full numpy arrays, so we can use allclose() and ignore
        # ordering of items with the same similarity value
        expected = matutils.sparse2full(expected, len(index))
        if num_best is not None: # when num_best is None, sims is already a numpy array
            sims = matutils.sparse2full(sims, len(index))
        self.assertTrue(numpy.allclose(expected, sims))
        if self.cls == similarities.Similarity:
            index.destroy()

        def test_shared_do_alias(self):
            dtype = self.dtype
            if dtype is None:
                dtype = theano.config.floatX

            rng = np.random.RandomState(utt.fetch_seed())
            x = np.asarray(rng.uniform(1, 2, [4, 2]), dtype=dtype)
            x = self.cast_value(x)
            x_ref = self.ref_fct(x)

            x_shared = self.shared_constructor(x, borrow=True)

            total = self.theano_fct(x_shared)

            total_func = theano.function([], total)

            total_val = total_func()

            assert np.allclose(self.ref_fct(x), total_val)

            x /= .5

            # not required by the contract but it is a feature we've implemented
            if self.shared_borrow_true_alias:
                assert np.allclose(self.ref_fct(x), total_func())
            else:
                assert np.allclose(x_ref, total_func())

def test_tnu():
    cosmo = core.FlatLambdaCDM(70.4, 0.272, Tcmb0=3.0)
    assert np.allclose(cosmo.Tnu0.value, 2.1412975665108247, rtol=1e-6)
    assert np.allclose(cosmo.Tnu(2).value, 6.423892699532474, rtol=1e-6)
    z = [0.0, 1.0, 2.0, 3.0]
    assert np.allclose(cosmo.Tnu(z), [2.14129757, 4.28259513,
                                      6.4238927, 8.56519027], rtol=1e-6)

 def testSymmetry(self):
     """Verify that the projection is symmetrical about the equator
     """
     for minDec in (-5.0, -1.0, 0.5):
         maxDec = minDec + 2.0
         config = EquatSkyMap.ConfigClass()
         config.decRange = minDec, maxDec
         skyMap = EquatSkyMap(config)
         for tractInfo in skyMap[0:1]:
             numPatches = tractInfo.getNumPatches()
             midXIndex = numPatches[0] / 2
             minPixelPosList = []
             maxPixelPosList = []
             maxYInd = numPatches[1] - 1
             for xInd in (0, midXIndex, numPatches[0] - 1):
                 minDecPatchInfo = tractInfo.getPatchInfo((xInd,0))
                 minDecPosBox = afwGeom.Box2D(minDecPatchInfo.getOuterBBox())
                 minPixelPosList += [
                     minDecPosBox.getMin(),
                     afwGeom.Point2D(minDecPosBox.getMaxX(), minDecPosBox.getMinY()),
                 ]
                 
                 maxDecPatchInfo = tractInfo.getPatchInfo((xInd, maxYInd))
                 maxDecPosBox = afwGeom.Box2D(maxDecPatchInfo.getOuterBBox())
                 maxPixelPosList += [
                     maxDecPosBox.getMax(),
                     afwGeom.Point2D(maxDecPosBox.getMinX(), maxDecPosBox.getMaxY()),
                 ]
             wcs = tractInfo.getWcs()
             minDecList = [wcs.pixelToSky(pos).getPosition(afwGeom.degrees)[1] for pos in minPixelPosList]
             maxDecList = [wcs.pixelToSky(pos).getPosition(afwGeom.degrees)[1] for pos in maxPixelPosList]
             self.assertTrue(numpy.allclose(minDecList, minDecList[0]))
             self.assertTrue(numpy.allclose(maxDecList, maxDecList[0]))
             self.assertTrue(minDecList[0] <= minDec)
             self.assertTrue(maxDecList[0] >= maxDec)

def test_warp_no_reproject_bounds_res(runner, tmpdir):
    srcname = 'tests/data/shade.tif'
    outputname = str(tmpdir.join('test.tif'))
    out_bounds = [-11850000, 4810000, -11849000, 4812000]
    result = runner.invoke(main_group, [
        'warp', srcname, outputname, '--res', 30, '--bounds'] + out_bounds)
    assert result.exit_code == 0
    assert os.path.exists(outputname)

    with rasterio.open(srcname) as src:
        with rasterio.open(outputname) as output:
            assert output.crs == src.crs
            assert np.allclose(output.bounds, out_bounds)
            assert np.allclose([30, 30], [output.transform.a, -output.transform.e])
            assert output.width == 34
            assert output.height == 67

    # dst-bounds should be an alias to bounds
    outputname = str(tmpdir.join('test2.tif'))
    out_bounds = [-11850000, 4810000, -11849000, 4812000]
    result = runner.invoke(main_group, [
        'warp', srcname, outputname, '--res', 30, '--dst-bounds'] + out_bounds)
    assert result.exit_code == 0
    assert os.path.exists(outputname)
    with rasterio.open(srcname) as src:
        with rasterio.open(outputname) as output:
            assert np.allclose(output.bounds, out_bounds)

def test_connected(Simulator):
    m = nengo.Network(label='test_connected', seed=123)
    with m:
        input = nengo.Node(output=lambda t: np.sin(t), label='input')
        output = nengo.Node(output=lambda t, x: np.square(x),
                            size_in=1,
                            label='output')
        nengo.Connection(input, output, synapse=None)  # Direct connection
        p_in = nengo.Probe(input, 'output')
        p_out = nengo.Probe(output, 'output')

    sim = Simulator(m)
    runtime = 0.5
    sim.run(runtime)

    with Plotter(Simulator) as plt:
        t = sim.trange()
        plt.plot(t, sim.data[p_in], label='sin')
        plt.plot(t, sim.data[p_out], label='sin squared')
        plt.plot(t, np.sin(t), label='ideal sin')
        plt.plot(t, np.sin(t) ** 2, label='ideal squared')
        plt.legend(loc='best')
        plt.savefig('test_node.test_connected.pdf')
        plt.close()

    sim_t = sim.trange()
    sim_sin = sim.data[p_in].ravel()
    sim_sq = sim.data[p_out].ravel()
    t = 0.001 * np.arange(len(sim_t))

    assert np.allclose(sim_t, t)
    # 1-step delay
    assert np.allclose(sim_sin[1:], np.sin(t[:-1]))
    assert np.allclose(sim_sq[1:], sim_sin[:-1] ** 2)

def test_normalization():
    """Test that `match_template` gives the correct normalization.

    Normalization gives 1 for a perfect match and -1 for an inverted-match.
    This test adds positive and negative squares to a zero-array and matches
    the array with a positive template.
    """
    n = 5
    N = 20
    ipos, jpos = (2, 3)
    ineg, jneg = (12, 11)
    image = np.full((N, N), 0.5)
    image[ipos:ipos + n, jpos:jpos + n] = 1
    image[ineg:ineg + n, jneg:jneg + n] = 0

    # white square with a black border
    template = np.zeros((n + 2, n + 2))
    template[1:1 + n, 1:1 + n] = 1

    result = match_template(image, template)

    # get the max and min results.
    sorted_result = np.argsort(result.flat)
    iflat_min = sorted_result[0]
    iflat_max = sorted_result[-1]
    min_result = np.unravel_index(iflat_min, result.shape)
    max_result = np.unravel_index(iflat_max, result.shape)

    # shift result by 1 because of template border
    assert np.all((np.array(min_result) + 1) == (ineg, jneg))
    assert np.all((np.array(max_result) + 1) == (ipos, jpos))

    assert np.allclose(result.flat[iflat_min], -1)
    assert np.allclose(result.flat[iflat_max], 1)

def isSkewSymmetric(A):
    """
        Returns True if input matrix is skew-symmetric.
        
        Parameter
        ---------
        A : array-like
            The input matrix. 
        
        Returns
        -------
        isSkewSymmetric : bool
            Returns True if the matrix is skew-symmetric; False otherwise.
        
        See Also
        --------
        isSquare, isSymmetric, isUpperTriangular, isLowerTriangular, isDiagonal
    """
    
    assert isSquare(A), 'Input matrix should be a square matrix.'
    
    if np.allclose(A, -1*A.T) and np.allclose(A.diagonal(),0):
        return True
    else:
        return False

def test_constant_scalar(Simulator, nl, plt, seed):
    """A Network that represents a constant value."""
    N = 30
    val = 0.5

    m = nengo.Network(label='test_constant_scalar', seed=seed)
    with m:
        m.config[nengo.Ensemble].neuron_type = nl()
        input = nengo.Node(output=val, label='input')
        A = nengo.Ensemble(N, 1)
        nengo.Connection(input, A)
        in_p = nengo.Probe(input, 'output')
        A_p = nengo.Probe(A, 'decoded_output', synapse=0.1)

    sim = Simulator(m, dt=0.001)
    sim.run(1.0)

    t = sim.trange()
    plt.plot(t, sim.data[in_p], label='Input')
    plt.plot(t, sim.data[A_p], label='Neuron approximation, pstc=0.1')
    plt.ylim([0, 1.05 * val])
    plt.legend(loc=0)

    assert np.allclose(sim.data[in_p], val, atol=.1, rtol=.01)
    assert np.allclose(sim.data[A_p][-10:], val, atol=.1, rtol=.01)

 def test_hv(self):
     def fn(x):
         return np.dot(x, x).sum()
     x = np.ones((3, 3))
     F = HessianVector(fn)
     self.assertTrue(np.allclose(x * 6, F(x, vectors=x)))
     self.assertTrue(np.allclose(x * 2, F(x[0], vectors=x[0])))

def _test_pes(Simulator, nl, plt, seed,
              pre_neurons=False, post_neurons=False, weight_solver=False,
              vin=np.array([0.5, -0.5]), vout=None, n=200,
              function=None, transform=np.array(1.), rate=1e-3):

    vout = np.array(vin) if vout is None else vout

    model = nengo.Network(seed=seed)
    with model:
        model.config[nengo.Ensemble].neuron_type = nl()

        u = nengo.Node(output=vin)
        v = nengo.Node(output=vout)
        a = nengo.Ensemble(n, dimensions=u.size_out)
        b = nengo.Ensemble(n, dimensions=u.size_out)
        e = nengo.Ensemble(n, dimensions=v.size_out)

        nengo.Connection(u, a)

        bslice = b[:v.size_out] if v.size_out < u.size_out else b
        pre = a.neurons if pre_neurons else a
        post = b.neurons if post_neurons else bslice

        conn = nengo.Connection(pre, post,
                                function=function, transform=transform,
                                learning_rule_type=PES(rate))
        if weight_solver:
            conn.solver = nengo.solvers.LstsqL2(weights=True)

        nengo.Connection(v, e, transform=-1)
        nengo.Connection(bslice, e)
        nengo.Connection(e, conn.learning_rule)

        b_p = nengo.Probe(bslice, synapse=0.03)
        e_p = nengo.Probe(e, synapse=0.03)

        weights_p = nengo.Probe(conn, 'weights', sample_every=0.01)
        corr_p = nengo.Probe(conn.learning_rule, 'correction', synapse=0.03)

    with Simulator(model) as sim:
        sim.run(0.5)
    t = sim.trange()
    weights = sim.data[weights_p]

    plt.subplot(311)
    plt.plot(t, sim.data[b_p])
    plt.ylabel("Post decoded value")
    plt.subplot(312)
    plt.plot(t, sim.data[e_p])
    plt.ylabel("Error decoded value")
    plt.subplot(313)
    plt.plot(t, sim.data[corr_p] / rate)
    plt.ylabel("PES correction")
    plt.xlabel("Time (s)")

    tend = t > 0.4
    assert np.allclose(sim.data[b_p][tend], vout, atol=0.05)
    assert np.allclose(sim.data[e_p][tend], 0, atol=0.05)
    assert np.allclose(sim.data[corr_p][tend] / rate, 0, atol=0.05)
    assert not np.allclose(weights[0], weights[-1], atol=1e-5)

def test_warp_reproject_dst_bounds(runner, tmpdir):
    """--bounds option works."""
    srcname = 'tests/data/shade.tif'
    outputname = str(tmpdir.join('test.tif'))
    out_bounds = [-106.45036, 39.6138, -106.44136, 39.6278]
    result = runner.invoke(
        main_group, [
            'warp', srcname, outputname, '--dst-crs', 'EPSG:4326',
            '--res', 0.001, '--bounds'] + out_bounds)
    assert result.exit_code == 0
    assert os.path.exists(outputname)

    with rasterio.open(outputname) as output:
        assert output.crs == {'init': 'epsg:4326'}
        assert np.allclose(output.bounds[0::3],
                           [-106.45036, 39.6278])
        assert np.allclose([0.001, 0.001],
                           [output.transform.a, -output.transform.e])

        # XXX: an extra row and column is produced in the dataset
        # because we're using ceil instead of floor internally.
        # Not necessarily a bug, but may change in the future.
        assert np.allclose([output.bounds[2] - 0.001, output.bounds[1] + 0.001],
                           [-106.44136, 39.6138])
        assert output.width == 10
        assert output.height == 15

        def test_return_internal_type(self):
            dtype = self.dtype
            if dtype is None:
                dtype = theano.config.floatX

            rng = np.random.RandomState(utt.fetch_seed())
            x = np.asarray(rng.uniform(0, 1, [2, 4]), dtype=dtype)
            x = self.cast_value(x)

            x_ref = self.ref_fct(x)
            x_shared = self.shared_constructor(x, borrow=False)
            total = self.theano_fct(x_shared)

            total_func = theano.function([], total)

            # in this case we can alias with the internal value
            x = x_shared.get_value(borrow=True, return_internal_type=True)
            assert self.test_internal_type(x)

            x /= .5

            # this is not required by the contract but it is a feature we can
            # implement for some type of SharedVariable.
            assert np.allclose(self.ref_fct(x), total_func())

            x = x_shared.get_value(borrow=False, return_internal_type=True)
            assert self.test_internal_type(x)
            assert x is not x_shared.container.value
            x /= .5

            # this is required by the contract
            assert not np.allclose(self.ref_fct(x), total_func())

def svd_example():
    a = np.floor(np.random.rand(4, 4)*20-6)
    logger.info("Matrix A:\n %s", a)
    b = np.floor(np.random.rand(4, 1)*20-6)
    logger.info("Matrix B:\n %s", b)

    u, s, v_t = np.linalg.svd(a) # SVD decomposition of A
    logger.info("Matrix U:\n %s", u)
    logger.info("Matrix S:\n %s", s)
    logger.info("Matrix V(transpose:\n %s", u)

    logger.info("Computing inverse using linalg.pinv")
    # Computing the inverse using pinv
    inv_pinv = np.linalg.pinv(a)
    logger.info("pinv:\n %s", inv_pinv)

    # Computing inverse using matrix decomposition
    logger.info("Computing inverse using svd matrix decomposition")
    inv_svd = np.dot(np.dot(v_t.T, np.linalg.inv(np.diag(s))), u.T)
    logger.info("svd inverse:\n %s", inv_svd)
    logger.info("comparing the results from pinv and svd_inverse:\n %s",
                np.allclose(inv_pinv, inv_svd))

    logger.info("Sol1: Solving x using pinv matrix... x=A^-1 x b")
    result_pinv_x = np.dot(inv_pinv, b)

    logger.info("Sol2: Solving x using svd_inverse matrix... x=A^-1 x b")
    result_svd_x = np.dot(inv_svd, b)

    if not np.allclose(result_pinv_x, result_svd_x):
        raise ValueError('Should have been True')

def test_tcmb():
    cosmo = core.FlatLambdaCDM(70.4, 0.272, Tcmb0=3.0)
    assert np.allclose(cosmo.Tcmb0.value, 3.0)
    assert np.allclose(cosmo.Tcmb(2).value, 9.0)
    z = [0.0, 1.0, 2.0, 3.0, 9.0]
    assert np.allclose(cosmo.Tcmb(z).value,
                       [3.0, 6.0, 9.0, 12.0, 30.0], rtol=1e-6)

    def test_bprop(self):
        r = []

        context = Context()
        for i in xrange(self.N):
            a = self.get_random_array()
            a_gpu = Connector(GpuMatrix.from_npa(a, 'float'), bu_device_id=context)
            vpooling_block = MeanPoolingBlock(a_gpu, axis=0)
            voutput, dL_dvoutput = vpooling_block.output.register_usage(context, context)
            _dL_voutput = self.get_random_array((dL_dvoutput.nrows, dL_dvoutput.ncols))
            GpuMatrix.from_npa(_dL_voutput, 'float').copy_to(context, dL_dvoutput)

            hpooling_block = MeanPoolingBlock(a_gpu, axis=1)
            houtput, dL_dhoutput = hpooling_block.output.register_usage(context, context)
            _dL_houtput = self.get_random_array((dL_dhoutput.nrows, dL_dhoutput.ncols))
            GpuMatrix.from_npa(_dL_houtput, 'float').copy_to(context, dL_dhoutput)

            vpooling_block.fprop()
            vpooling_block.bprop()
            dL_dmatrix = vpooling_block.dL_dmatrix.to_host()
            r.append(np.allclose(dL_dmatrix,
                                 np.repeat(_dL_voutput/a.shape[0], a.shape[0], 0),
                                 atol=1e-6))

            hpooling_block.fprop()
            hpooling_block.bprop()
            hpooling_block.dL_dmatrix.to_host()
            dL_dmatrix = hpooling_block.dL_dmatrix.to_host()
            r.append(np.allclose(dL_dmatrix,
                                 np.repeat(_dL_houtput/a.shape[1], a.shape[1], 1),
                                 atol=1e-6))

        self.assertEqual(sum(r), 2 * self.N)