def _testUniformSampleMultiDimensional(self):
    # DISABLED: Please enable this test once b/issues/30149644 is resolved.
    with self.test_session():
      batch_size = 2
      a_v = [3.0, 22.0]
      b_v = [13.0, 35.0]
      a = constant_op.constant([a_v] * batch_size)
      b = constant_op.constant([b_v] * batch_size)

      uniform = uniform_lib.Uniform(low=a, high=b)

      n_v = 100000
      n = constant_op.constant(n_v)
      samples = uniform.sample(n)
      self.assertEqual(samples.get_shape(), (n_v, batch_size, 2))

      sample_values = self.evaluate(samples)

      self.assertFalse(
          np.any(sample_values[:, 0, 0] < a_v[0]) or
          np.any(sample_values[:, 0, 0] >= b_v[0]))
      self.assertFalse(
          np.any(sample_values[:, 0, 1] < a_v[1]) or
          np.any(sample_values[:, 0, 1] >= b_v[1]))

      self.assertAllClose(
          sample_values[:, 0, 0].mean(), (a_v[0] + b_v[0]) / 2, atol=1e-2)
      self.assertAllClose(
          sample_values[:, 0, 1].mean(), (a_v[1] + b_v[1]) / 2, atol=1e-2)

 def testDtype(self):
   with self.test_session():
     d = array_ops.fill([2, 3], 12., name="fill")
     self.assertEqual(d.get_shape(), [2, 3])
     # Test default type for both constant size and dynamic size
     z = array_ops.zeros([2, 3])
     self.assertEqual(z.dtype, dtypes_lib.float32)
     self.assertEqual([2, 3], z.get_shape())
     self.assertAllEqual(z.eval(), np.zeros([2, 3]))
     z = array_ops.zeros(array_ops.shape(d))
     self.assertEqual(z.dtype, dtypes_lib.float32)
     self.assertEqual([2, 3], z.get_shape())
     self.assertAllEqual(z.eval(), np.zeros([2, 3]))
     # Test explicit type control
     for dtype in [
         dtypes_lib.float32, dtypes_lib.float64, dtypes_lib.int32,
         dtypes_lib.uint8, dtypes_lib.int16, dtypes_lib.int8,
         dtypes_lib.complex64, dtypes_lib.complex128, dtypes_lib.int64,
         dtypes_lib.bool, dtypes_lib.string
     ]:
       z = array_ops.zeros([2, 3], dtype=dtype)
       self.assertEqual(z.dtype, dtype)
       self.assertEqual([2, 3], z.get_shape())
       z_value = z.eval()
       self.assertFalse(np.any(z_value))
       self.assertEqual((2, 3), z_value.shape)
       z = array_ops.zeros(array_ops.shape(d), dtype=dtype)
       self.assertEqual(z.dtype, dtype)
       self.assertEqual([2, 3], z.get_shape())
       z_value = z.eval()
       self.assertFalse(np.any(z_value))
       self.assertEqual((2, 3), z_value.shape)

    def tag_sites(self, scaled_positions, symprec=1e-3):
        """Returns an integer array of the same length as *scaled_positions*, 
        tagging all equivalent atoms with the same index.

        Example:

        >>> from ase.lattice.spacegroup import Spacegroup
        >>> sg = Spacegroup(225)  # fcc
        >>> sg.tag_sites([[0.0, 0.0, 0.0], 
        ...               [0.5, 0.5, 0.0], 
        ...               [1.0, 0.0, 0.0], 
        ...               [0.5, 0.0, 0.0]])
        array([0, 0, 0, 1])
        """
        scaled = np.array(scaled_positions, ndmin=2)
        scaled %= 1.0
        scaled %= 1.0
        tags = -np.ones((len(scaled), ), dtype=int)
        mask = np.ones((len(scaled), ), dtype=np.bool)
        rot, trans = self.get_op()
        i = 0
        while mask.any():
            pos = scaled[mask][0]
            sympos = np.dot(rot, pos) + trans
            # Must be done twice, see the scaled_positions.py test
            sympos %= 1.0
            sympos %= 1.0
            m = ~np.all(np.any(np.abs(scaled[np.newaxis,:,:] - 
                                      sympos[:,np.newaxis,:]) > symprec, 
                               axis=2), axis=0)
            assert not np.any((~mask) & m)
            tags[m] = i
            mask &= ~m
            i += 1
        return tags

def allowed_region( V_nj, ave_j ):

    # read PCs
    PC1 = V_nj[0]
    PC2 = V_nj[1]
    n_band = len( PC1 )
    band_ticks = np.arange( n_band )

    x_ticks = np.linspace(-0.4,0.2,RESOLUTION)
    y_ticks = np.linspace(-0.2,0.4,RESOLUTION)
    x_mesh, y_mesh, band_mesh = np.meshgrid( x_ticks, y_ticks, band_ticks, indexing='ij' )
    vec_mesh = x_mesh * PC1[ band_mesh ] + y_mesh * PC2[ band_mesh ] + ave_j[ band_mesh ]

    x_grid, y_grid = np.meshgrid( x_ticks, y_ticks, indexing='ij' )
    prohibited_grid = np.zeros_like( x_grid )

    for ii in xrange( len( x_ticks ) ) :
        for jj in xrange( len( y_ticks ) ) :

            if np.any( vec_mesh[ii][jj] < 0. ) :
                prohibited_grid[ii][jj] = 1
                if np.any( vec_mesh[ii][jj] > 1. ) :
                    prohibited_grid[ii][jj] = 3
            elif np.any( vec_mesh[ii][jj] > 1. ) :
                prohibited_grid[ii][jj] = 2
            else :
                prohibited_grid[ii][jj] = 0

    return x_grid, y_grid, prohibited_grid

 def predict(self, session, X, y=None):
   """Make predictions from the provided model."""
   # If y is given, the loss is also calculated
   # We deactivate dropout by setting it to 1
   dp = 1
   losses = []
   results = []
   if np.any(y):
       data = data_iterator(X, y, batch_size=self.config.batch_size,
                            label_size=self.config.label_size, shuffle=False)
   else:
       data = data_iterator(X, batch_size=self.config.batch_size,
                            label_size=self.config.label_size, shuffle=False)
   for step, (x, y) in enumerate(data):
     feed = self.create_feed_dict(input_batch=x, dropout=dp)
     if np.any(y):
       feed[self.labels_placeholder] = y
       loss, preds = session.run(
           [self.loss, self.predictions], feed_dict=feed)
       losses.append(loss)
     else:
       preds = session.run(self.predictions, feed_dict=feed)
     predicted_indices = preds.argmax(axis=1)
     results.extend(predicted_indices)
   return np.mean(losses), results

    def test_using_gpu_3(self):

        if theano.config.device.find('gpu') > -1:

            from theano import function, config, shared, sandbox, Out
            import theano.tensor as T
            import numpy
            import time

            vlen = 10 * 30 * 70  # 10 x #cores x # threads per core
            iters = 10

            rng = numpy.random.RandomState(22)
            x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
            f = function([],
                    Out(sandbox.cuda.basic_ops.gpu_from_host(T.exp(x)),
                        borrow=True))
            # print f.maker.fgraph.toposort()
            t0 = time.time()
            for i in xrange(iters):
                r = f()
            t1 = time.time()
            print 'Looping %d times took' % iters, t1 - t0, 'seconds'
            print 'Result is', r
            print 'Numpy result is', numpy.asarray(r)
            if numpy.any([isinstance(x.op, T.Elemwise)
                          for x in f.maker.fgraph.toposort()]):
                print 'Used the cpu'
            else:
                print 'Used the gpu'

            assert not numpy.any([isinstance(x.op, T.Elemwise)
                                  for x in f.maker.fgraph.toposort()])

def fill_betweenx_discontinuous(ax, ymin, ymax, x, freq=1, **kwargs):
    """Fill betwwen x even if x is discontinuous clusters
    Parameters
    ----------
    ax : axis
    x : list

    Returns
    -------
    ax : axis
    """
    x = np.array(x)
    min_gap = 1.1 / freq
    while np.any(x):
        # If with single time point
        if len(x) > 1:
            xmax = np.where((x[1:] - x[:-1]) > min_gap)[0]
        else:
            xmax = [0]

        # If continuous
        if not np.any(xmax):
            xmax = [len(x) - 1]

        ax.fill_betweenx((ymin, ymax), x[0], x[xmax[0]], **kwargs)

        # remove from list
        x = x[(xmax[0] + 1) :]
    return ax

def get_polar_motion(time):
    """
    gets the two polar motion components in radians for use with apio13
    """
    # Get the polar motion from the IERS table
    xp, yp, status = iers.IERS_Auto.open().pm_xy(time, return_status=True)

    wmsg = None
    if np.any(status == iers.TIME_BEFORE_IERS_RANGE):
        wmsg = ('Tried to get polar motions for times before IERS data is '
                'valid. Defaulting to polar motion from the 50-yr mean for those. '
                'This may affect precision at the 10s of arcsec level')
        xp.ravel()[status.ravel() == iers.TIME_BEFORE_IERS_RANGE] = _DEFAULT_PM[0]
        yp.ravel()[status.ravel() == iers.TIME_BEFORE_IERS_RANGE] = _DEFAULT_PM[1]

        warnings.warn(wmsg, AstropyWarning)

    if np.any(status == iers.TIME_BEYOND_IERS_RANGE):
        wmsg = ('Tried to get polar motions for times after IERS data is '
                'valid. Defaulting to polar motion from the 50-yr mean for those. '
                'This may affect precision at the 10s of arcsec level')

        xp.ravel()[status.ravel() == iers.TIME_BEYOND_IERS_RANGE] = _DEFAULT_PM[0]
        yp.ravel()[status.ravel() == iers.TIME_BEYOND_IERS_RANGE] = _DEFAULT_PM[1]

        warnings.warn(wmsg, AstropyWarning)

    return xp.to(u.radian).value, yp.to(u.radian).value

 def backward(self, top, propagate_down, bottom):
     h=bottom[0].data.shape[2]
     w=bottom[0].data.shape[3]
     num_proposals = bottom[1].data.shape[0]
     fbox = self.fbox
     bottom[0].diff[...] = 0
     for b in np.arange(num_proposals):
         if np.any(self.act[b]==True):
             wb=fbox[b,2]-fbox[b,0]
             hb=fbox[b,3]-fbox[b,1]
             if hb==0 or wb==0: #bounding box with size 0
                 pass
             elif h==hb and w==wb: #bounding box with size of the image
                 pass
             else:
                 bottom[0].diff[0,self.act[b],:,:] += bottom[1].data[b,self.act[b],np.newaxis,np.newaxis]/(h*w-hb*wb) #negative
                 bottom[0].diff[0,self.act[b],fbox[b,1]:fbox[b,3]+1,fbox[b,0]:fbox[b,2]+1] += -bottom[1].data[b,self.act[b],np.newaxis,np.newaxis]/(h*w-hb*wb) #negative
                 bottom[0].diff[0,self.act[b],fbox[b,1]:fbox[b,3]+1,fbox[b,0]:fbox[b,2]+1] += -bottom[1].data[b,self.act[b],np.newaxis,np.newaxis]/(hb*wb) #positive
     bottom[0].diff[...] *= self.myloss_weight
     bottom[1].diff[...] = self.myloss_weight*self.hinge
     if np.any(np.isnan(bottom[0].diff)):
         print "Nan error"
         dsfsf
     if np.any(np.isnan(bottom[1].diff)):
         print "Nan error"
         dsfsf

    def sample_representer_points(self):
        # Sample representer points only in the
        # configuration space by setting all environmental
        # variables to 1
        D = np.where(self.is_env == 0)[0].shape[0]

        lower = self.lower[np.where(self.is_env == 0)]
        upper = self.upper[np.where(self.is_env == 0)]

        self.sampling_acquisition.update(self.model)

        for i in range(5):
            restarts = np.random.uniform(low=lower,
                                         high=upper,
                                         size=(self.Nb, D))
            sampler = emcee.EnsembleSampler(self.Nb, D,
                                        self.sampling_acquisition_wrapper)

            self.zb, self.lmb, _ = sampler.run_mcmc(restarts, 50)
            if not np.any(np.isinf(self.lmb)):
                break
            else:
                print("Infinity")
        if np.any(np.isinf(self.lmb)):
            raise ValueError("Could not sample valid representer points! LogEI is -infinity")
        if len(self.zb.shape) == 1:
            self.zb = self.zb[:, None]
        if len(self.lmb.shape) == 1:
            self.lmb = self.lmb[:, None]

        # Project representer points to subspace
        proj = np.ones([self.zb.shape[0],
                    self.upper[self.is_env == 1].shape[0]])
        proj *= self.upper[self.is_env == 1].shape[0]
        self.zb = np.concatenate((self.zb, proj), axis=1)

def algebraic2parametric(coeff):
        '''
        Based on matlab function "ellipse_param.m" which accompanies
        "Least-Squares Fitting of Circles and Ellipses", W. Gander, G. H. Golub, R. Strebel,
                BIT Numerical Mathematics, Springer 1994

        convert the coefficients (a,b,c,d,e,f) of the algebraic equation:
                ax^2 + bxy + cy^2 + dx + ey + f = 0
        to the parameters of the parametric equation.  The parameters are
        returned as a dictionary containing:
                center - center of the ellipse
                a - major axis
                b - minor axis
                alpha - angle of major axis

        convention note: alpha is measured as positive values towards the y-axis
        '''
        #print coeff
        #print ("A=%.3f B=%.3f C=%.3f D=%.3f E=%.3f F=%.3f"
        #       %(coeff[0], coeff[1], coeff[2], coeff[3], coeff[4], coeff[5],))

        if numpy.any(numpy.isnan(coeff)) or numpy.any(numpy.isinf(coeff)):
                return None

        A   = numpy.array((coeff[0], coeff[1]/2, coeff[1]/2, coeff[2]))
        A.shape = 2,2
        bb  = numpy.asarray(coeff[3:5])
        c   = coeff[5]

        D,Q = scipy.linalg.eig(A)
        D = D.real
        det = D[0]*D[1]
        if det <= 0:
                return None
        else: 
                bs = numpy.dot(Q.transpose(), bb)
                alpha = numpy.arctan2(Q[1,0], Q[0,0])

                zs = scipy.linalg.solve(-2*numpy.diagflat(D), bs)
                z = numpy.dot(Q, zs)
                h = numpy.dot(-bs.transpose(), zs) / 2 - c

                a = numpy.sqrt(h/D[0])
                b = numpy.sqrt(h/D[1])

        ## correct backwards major/minor axes
        ## 'major axis as a, minor axis as b'
        if b > a:
                temp = b
                b = a
                a = temp
                alpha = math.pi/2 + alpha

        #print "alpha", alpha
        if alpha <= -math.pi/2:
                alpha += math.pi
        elif alpha > math.pi/2:
                alpha -= math.pi

        return {'center':z, 'a':a, 'b':b, 'alpha':alpha}

    def collide(self):
        for a in self.actors:
            a.collision_prepare()

        for i, ai in enumerate(self.actors):
            if isinstance(ai, ParticleActor):
                ai.collision_self()

            for j, aj in enumerate(self.actors):
                if isinstance(aj, ParticleActor):
                    continue
                if i == j: continue
                info = spatial.Info(ai.spatial_grid, aj.spatial_mesh, i==j)
                mask = info.triangle != -1
                active = np.flatnonzero(mask)   # active vertex idx
                if np.any(mask):# and not isinstance(ai, StaticActor):
                    triangle = info.triangle[active]
                    bary = info.bary[active]

                    velocity = aj.velocity[aj.mesh.faces[triangle]]
                    velocity = np.einsum('vtc,vt->vc', velocity, bary)
                    relative_velocity = ai.velocity[active] - velocity

                    friction = 1e-2
                    stiffness = 1e1
                    force = info.depth[active][:, None] * info.normal[active] * stiffness - relative_velocity * friction
                    assert not np.any(np.isnan(force))

                    np.add.at(ai.force, active, force)
                    corners = aj.mesh.faces[triangle]
                    for i in range(3):
                        np.add.at(aj.force, corners[:, i], -bary[:, [i]] * force)

    def testBinds(self):
        ds = normalFeatureDataset()
        ds_data = ds.samples.copy()
        ds_chunks = ds.chunks.copy()
        self.failUnless(N.all(ds.samples == ds_data)) # sanity check

        funcs = ['zscore', 'coarsenChunks']
        if externals.exists('scipy'):
            funcs.append('detrend')

        for f in funcs:
            eval('ds.%s()' % f)
            self.failUnless(N.any(ds.samples != ds_data) or
                            N.any(ds.chunks != ds_chunks),
                msg="We should have modified original dataset with %s" % f)
            ds.samples = ds_data.copy()
            ds.chunks = ds_chunks.copy()

        # and some which should just return results
        for f in ['aggregateFeatures', 'removeInvariantFeatures',
                  'getSamplesPerChunkLabel']:
            res = eval('ds.%s()' % f)
            self.failUnless(res is not None,
                msg='We should have got result from function %s' % f)
            self.failUnless(N.all(ds.samples == ds_data),
                msg="Function %s should have not modified original dataset" % f)

	def check_obs_scheme(self):
		" Checks the internal validity of provided observation schemes "

		# check sub_pops
		idx_union = np.sort(self._sub_pops[0])
		i = 1
		while idx_union.size < self._p and i < len(self._sub_pops):
			idx_union = np.union1d(idx_union, self._sub_pops[i]) 
			i += 1
		if idx_union.size != self._p or np.any(idx_union!=np.arange(self._p)):
			raise Exception(('all subpopulations together have to cover '
			'exactly all included observed varibles y_i in y.'
			'This is not the case. Change the difinition of '
			'subpopulations in variable sub_pops or reduce '
			'the number of observed variables p. '
			'The union of indices of all subpopulations is'),
			idx_union )

		# check obs_time
		if not self._obs_time[-1]==self._T:
			raise Exception(('Entries of obs_time give the respective ends of '
							'the periods of observation for any '
							'subpopulation. Hence the last entry of obs_time '
							'has to be the full recording length. The last '
							'entry of obs_time before is '), self._obs_time[-1])

		if np.any(np.diff(self._obs_time)<1):
			raise Exception(('lengths of observation have to be at least 1. '
							'Minimal observation time for a subpopulation: '),
							np.min(np.diff(self._obs_time)))

		# check obs_pops
		if not self._obs_time.size == self._obs_pops.size:
			raise Exception(('each entry of obs_pops gives the index of the '
							'subpopulation observed up to the respective '
							'time given in obs_time. Thus the sizes of the '
							'two arrays have to match. They do not. '
							'no. of subpop. switch points and no. of '
							'subpopulations ovserved up to switch points '
							'are '), (self._obs_time.size, self._obs_pops.size))

		idx_pops = np.sort(np.unique(self._obs_pops))
		if not np.min(idx_pops)==0:
			raise Exception(('first subpopulation has to have index 0, but '
							'is given the index '), np.min(idx_pops))
		elif not idx_pops.size == len(self._sub_pops):
			raise Exception(('number of specified subpopulations in variable '
							'sub_pops does not meet the number of '
							'subpopulations indexed in variable obs_pops. '
							'Delete subpopulations that are never observed, '
							'or change the observed subpopulations in '
							'variable obs_pops accordingly. The number of '
							'indexed subpopulations is '),
							len(self._sub_pops))
		elif not np.all(np.diff(idx_pops)==1):
			raise Exception(('subpopulation indices have to be consecutive '
							'integers from 0 to the total number of '
							'subpopulations. This is not the case. '
							'Given subpopulation indices are '),
							idx_pops)

    def pop_planes(geometry, kwargs):
        # Convert miller index specifications to normal vectors
        miller_defs = kwargs.pop("planes_miller", None)
        if miller_defs is not None:
            if np.any(np.all(abs(miller_defs[:,0:3]) < EPSILON, axis=1)):
                error("Emtpy miller index tuple")
            miller_defs[:,0:3] = miller_to_normal(
                np.dot(geometry.latvecs, geometry.bravais_cell),
                miller_defs[:,0:3])
        else:
            miller_defs = np.zeros((0, 4), dtype=float)
            
        # Convert plane normal vector specifications into cartesian coords.
        normal_defs = kwargs.pop("planes_normal", None)
        if normal_defs is not None:
            normal_defs[:,0:3] = geometry.coord_transform(
                normal_defs[:,0:3],
                kwargs.pop("planes_normal_coordsys", "lattice"))
            if np.any(np.all(abs(normal_defs[:,0:3]) < EPSILON, axis=1)):
                error("Emtpy normal vector definition")
        else:
            normal_defs = np.zeros((0, 4), dtype=float)

        # Append two defintions
        planes_normal = np.vstack(( miller_defs, normal_defs ))
        return planes_normal

def j_roots(n, alpha, beta, mu=0):
    """[x,w] = j_roots(n,alpha,beta)

    Returns the roots (x) of the nth order Jacobi polynomial, P^(alpha,beta)_n(x)
    and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting
    function (1-x)**alpha (1+x)**beta with alpha,beta > -1.
    """
    if any(alpha <= -1) or any(beta <= -1):
        raise ValueError("alpha and beta must be greater than -1.")
    assert n > 0, "n must be positive."

    (p, q) = (alpha, beta)
    # from recurrence relations
    sbn_J = (
        lambda k: 2.0
        / (2.0 * k + p + q)
        * sqrt((k + p) * (k + q) / (2 * k + q + p + 1))
        * (np.where(k == 1, 1.0, sqrt(k * (k + p + q) / (2.0 * k + p + q - 1))))
    )
    if any(p == q):  # XXX any or all???
        an_J = lambda k: 0.0 * k
    else:
        an_J = lambda k: np.where(
            k == 0, (q - p) / (p + q + 2.0), (q * q - p * p) / ((2.0 * k + p + q) * (2.0 * k + p + q + 2))
        )
    g = cephes.gamma
    mu0 = 2.0 ** (p + q + 1) * g(p + 1) * g(q + 1) / (g(p + q + 2))
    val = gen_roots_and_weights(n, an_J, sbn_J, mu0)
    if mu:
        return val + [mu0]
    else:
        return val

    def test_using_gpu_1(self):
        # I'm checking if this compiles and runs
        from theano import function, config, shared, sandbox
        import theano.tensor as T
        import numpy
        import time

        vlen = 10 * 30 * 70  # 10 x #cores x # threads per core
        iters = 10

        rng = numpy.random.RandomState(22)
        x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
        f = function([], T.exp(x))
        # print f.maker.fgraph.toposort()
        t0 = time.time()
        for i in xrange(iters):
            r = f()
        t1 = time.time()
        print 'Looping %d times took' % iters, t1 - t0, 'seconds'
        print 'Result is', r
        if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
            print 'Used the cpu'
        else:
            print 'Used the gpu'
        if theano.config.device.find('gpu') > -1:
            assert not numpy.any( [isinstance(x.op,T.Elemwise) for x in f.maker.fgraph.toposort()])
        else:
            assert numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()])

    def test_unwrap(self):
        """Test different geometry types are appropriately unwrapped."""

        wrapper = Wrapper()
        path = tempfile.mkdtemp()
        for desc, geom in self.possible.iteritems():
            unwrapped = wrapper.unwrap(geom)
            if desc in self.actual_unwrapped:
                self.assertTrue(self.actual_unwrapped[desc].almost_equals(unwrapped, decimal=5))
            try:
                self.assertEqual(type(geom), type(unwrapped))
            except AssertionError:
                if desc == 'axis_polygon':
                    # by necessity of being split on the axis, this will come out as a multipolygon
                    self.assertIsInstance(unwrapped, MultiPolygon)
                else:
                    raise

            self.assertFalse(np.any(np.array(unwrapped) < 0.0))
            if isinstance(unwrapped, (MultiPolygon, Polygon)):
                it = get_iter(unwrapped)
                for polygon in it:
                    self.assertFalse(np.any(np.array(polygon.exterior) > 360.0))
            else:
                self.assertFalse(np.any(np.array(unwrapped) > 360.0))

 def __init__(self, x, y):
             
     assert np.ndim(x)==2 and np.ndim(y)==2 and np.shape(x)==np.shape(y), \
         'x and y must be 2D arrays of the same size.'
     
     if np.any(np.isnan(x)) or np.any(np.isnan(y)):
         x = np.ma.masked_where( (isnan(x)) | (isnan(y)) , x)
         y = np.ma.masked_where( (isnan(x)) | (isnan(y)) , y)
         
     self.x_vert = x
     self.y_vert = y
     
     mask_shape = tuple([n-1 for n in self.x_vert.shape])
     self.mask_rho = np.ones(mask_shape, dtype='d')
     
     # If maskedarray is given for verticies, modify the mask such that 
     # non-existant grid points are masked.  A cell requires all four
     # verticies to be defined as a water point.
     if isinstance(self.x_vert, np.ma.MaskedArray):
         mask = (self.x_vert.mask[:-1,:-1] | self.x_vert.mask[1:,:-1] | \
                 self.x_vert.mask[:-1,1:] | self.x_vert.mask[1:,1:])
         self.mask_rho = np.asarray(~(~np.bool_(self.mask_rho) | mask), dtype='d')
     
     if isinstance(self.y_vert, np.ma.MaskedArray):
         mask = (self.y_vert.mask[:-1,:-1] | self.y_vert.mask[1:,:-1] | \
                 self.y_vert.mask[:-1,1:] | self.y_vert.mask[1:,1:])
         self.mask_rho = np.asarray(~(~np.bool_(self.mask_rho) | mask), dtype='d')
     
     self._calculate_subgrids()
     self._calculate_metrics()        

    def _LoadHeader(self):
        ConfigLoader._LoadHeader(self)
        bc = self.Domain.BlockCounts
        nBlocks = self.Domain.TotalBlocks
        # Number of blocks in its neighbourhood.
        self.BlockNeighbourhoodSize = np.zeros(nBlocks, dtype=np.uint8)
        # Number of blocks in the neighbourhood that are available
        self.BlockNeighbourhoodAvailable = np.zeros(nBlocks, dtype=np.uint8)
        # Number of blocks in the neighbourhood that are done
        self.BlockNeighbourhoodDone = np.zeros(nBlocks, dtype=np.uint8)
        # Is the block itself done
        self.IsBlockDone = np.zeros(nBlocks, dtype=np.bool)
        # Lock to ensure only one thread at a time updates IsBlockDone
        # and BlockNeighbourhoodDone
        self.DoneLock = threading.RLock()

        # Compute the size of each block's neighbourhood
        for bIjk, bIdx in self.Domain.BlockIndexer.IterBoth():
            for i, delta in enumerate(self.NeighbourhoodOffsets):
                nIdx = bIdx + delta
                if np.any(nIdx < 0) or np.any(nIdx >= bc):
                    continue
                
                nIjk = self.Domain.BlockIndexer.NdToOne(nIdx)
                self.BlockNeighbourhoodSize[nIjk] += 1
                continue
            continue
        return