def assertArraysEqual(self, arr1, arr2):
        """
        Ensure that two numpy / numpy.ma arrays are equivalent in both their
        mask and their data.
        """
        if (arr1.shape != arr2.shape):
            msg = "Shapes differ:\n" + \
              str(arr1.shape) + " != " + str(arr2.shape)
            raise AssertionError(msg)

        mask1 = ma.getmaskarray(arr1)
        mask2 = ma.getmaskarray(arr2)
        masks_equal = np.array_equal(mask1, mask2)
        if (not masks_equal):
            msg = "Masks differ:\n" + \
              str(mask1) + " != " + str(mask2) + "\n" + \
              "Arrays are: \n" + \
              str(arr1) + "\n" + \
              str(arr2)

            raise AssertionError(msg)

        vals_equal = ma.allequal(arr1, arr2)
        if (not vals_equal):
            msg = "Values differ:\n" + \
              str(arr1) + " != " + str(arr2)
            raise AssertionError(msg)

    def _sanitise_geometry(self, lons, lats):
        """
        Sanitise geometry by removing any masked points.
        :param lons:
        :param lats:
        :return: A tuple containing arrays of (lon, lat)
        """
        # Align the arrays
        lons, lats = self.__align_lons_lats(lons, lats)

        # Get masks
        lon_mask = ma.getmaskarray(lons)
        lat_mask = ma.getmaskarray(lats)

        # Filter the arrays
        # kltsa 19/07/2016 Change for issue 23330: Some filters 
        #                  were excluded in order to allow values 
        #                  up to 90 to be included.
        self.longitudes = lons[
            (lons >= -180) &
            (lons <= 180) &
            #(lats >= 90) &
            #(lats <= 90) &
            (lon_mask == False)# &
            #(lat_mask == False)
        ]

        self.latitudes = lats[
            #(lons >= -180) &
            #(lons <= 180) &
            (lats >= -90) &
            (lats <= 90) &
            #(lon_mask == False) &
            (lat_mask == False)
        ]

    def _sanitise_geometry(self, lons, lats):
        """
        Sanitise geometry by removing any masked points.
        :param lons:
        :param lats:
        :return: A tuple containing arrays of (lon, lat)
        """
        # Align the arrays
        lons, lats = self.__align_lons_lats(lons, lats)

        # Get masks
        lon_mask = ma.getmaskarray(lons)
        lat_mask = ma.getmaskarray(lats)

        # Filter the arrays
        self.longitudes = lons[
            (lons >= -180) &
            (lons <= 180) &
            (lats >= -90) &
            (lats <= 90) &
            (lon_mask == False) &
            (lat_mask == False)
        ]

        self.latitudes = lats[
            (lons >= -180) &
            (lons <= 180) &
            (lats >= -90) &
            (lats <= 90) &
            (lon_mask == False) &
            (lat_mask == False)
        ]

def _tabulate_time_series(a):
    """
    Private function called by tabulate for flexible-dtype TimeSeries.
    """
    basedtype = a.dtype
    basenames = basedtype.names
    if basenames is None:
        _varshape = a._varshape
        if _varshape != ():
            pseudodtype = [('_dates', int),
                           ('_data',(basedtype, _varshape)),
                           ('_mask',(bool,_varshape))]
        else:
            pseudodtype = [('_dates', int),
                           ('_data', basedtype),
                           ('_mask', bool)]
        pseudo = itertools.izip(a._dates, a.filled(), ma.getmaskarray(a),)
    else:
        pseudodtype = [('_dates', int)]
        pseudodtype.extend([(fname,[('_data',ftype), ('_mask',bool)])
                            for (fname,ftype) in basedtype.descr])
        fields = [a[f] for f in basenames]
        pseudo = itertools.izip(a._dates,
                                *[zip(f.filled().flat, ma.getmaskarray(f).flat)
                                  for f in fields])
    return np.fromiter(pseudo, dtype=pseudodtype)

    def _assertMaskedArray(self, assertion, a, b, strict, **kwargs):
        # Define helper function to extract unmasked values as a 1d
        # array.
        def unmasked_data_as_1d_array(array):
            array = ma.asarray(array)
            if array.ndim == 0:
                if array.mask:
                    data = np.array([])
                else:
                    data = np.array([array.data])
            else:
                data = array.data[~ma.getmaskarray(array)]
            return data

        # Compare masks. This will also check that the array shapes
        # match, which is not tested when comparing unmasked values if
        # strict is False.
        a_mask, b_mask = ma.getmaskarray(a), ma.getmaskarray(b)
        np.testing.assert_array_equal(a_mask, b_mask)

        if strict:
            assertion(a.data, b.data, **kwargs)
        else:
            assertion(unmasked_data_as_1d_array(a),
                      unmasked_data_as_1d_array(b),
                      **kwargs)

def errorbar(ax, x, y, xerr=True, yerr=True, fmt='-', ecolor=None, elinewidth=None, capsize=3, 
             barsabove=False, lolims=False, uplims=False, xlolims=False, xuplims=False, **kwargs):
    """Plots two Dvects against each other. Requires matplotlib.
    
    See matplotlib.axes.Axes.errorbar for an explanation of most fields. The only differences are:

    ax   -- the Axes you are plotting in, needed to allow a call of the form ax.errorbar(...) to the
            standard matplotlib errorbar plotter
    x    -- x Dvect
    y    -- y Dvect

    Quite often you may want to set fmt=.' to get points. The lines are for compatibility with
    the matplotlib errorbar routine.
    """
    import warnings
    
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        raise DvectError("dvect.errorbar: matplotlib needed for plotting") 

    # Identify joint OK part
    ok = ~(getmaskarray(x) | getmaskarray(y))

    # Catch xerr and yerr which are re-defined cf maplotlib errorbar since 
    # Dvects carry their own errors
    xerr = x.err.data[ok] if xerr and x.err is not None and isinstance(x,Dvect) else None
    yerr = y.err.data[ok] if yerr and y.err is not None and isinstance(y,Dvect) else None

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        ax.errorbar(x.data[ok], y.data[ok], xerr=xerr, yerr=yerr, fmt=fmt, ecolor=ecolor, elinewidth=elinewidth, capsize=capsize, 
                    barsabove=barsabove, lolims=lolims, uplims=uplims, xlolims=xlolims, xuplims=xuplims, **kwargs)

 def _compute_masks_differ(self, var1, var2):
     if (np.array_equal(
         ma.getmaskarray(var1),
         ma.getmaskarray(var2))):
         return False
     else:
         return True

 def __setslice__(self, i, j, value):
     "Sets the slice described by [i,j] to `value`."
     _localdict = self.__dict__
     d = self._data
     m = _localdict['_fieldmask']
     names = self.dtype.names
     if value is masked:
         for n in names:
             m[i:j][n] = True
     elif not self._hardmask:
         fval = filled(value)
         mval = getmaskarray(value)
         for n in names:
             d[n][i:j] = fval
             m[n][i:j] = mval
     else:
         mindx = getmaskarray(self)[i:j]
         dval = np.asarray(value)
         valmask = getmask(value)
         if valmask is nomask:
             for n in names:
                 mval = mask_or(m[n][i:j], valmask)
                 d[n][i:j][~mval] = value
         elif valmask.size > 1:
             for n in names:
                 mval = mask_or(m[n][i:j], valmask)
                 d[n][i:j][~mval] = dval[~mval]
                 m[n][i:j] = mask_or(m[n][i:j], mval)
     self._fieldmask = m

def self_training(X, y, X_unLabeled, clf, th):
    clf.fit(X=X, y=y)
    index_unlabeled = ma.arange(0, len(X_unLabeled), 1)
    y_unlabeled = np.zeros(len(X_unLabeled))
    train_is_failed = False

    while True:
        probs = clf.predict_proba(X=X_unLabeled[~ma.getmaskarray(index_unlabeled)])
        index_greater_equal = np.greater_equal([max(d) for d in probs], [th]*len(probs))
        index_labelable = index_unlabeled.data[~ma.getmaskarray(index_unlabeled)][index_greater_equal]

        if not len(index_labelable) > 0:
            if not len(index_unlabeled.data[ma.getmaskarray(index_unlabeled)]) > 0:
                train_is_failed = True
            break

        index_unlabeled[index_labelable] = ma.masked

        if index_unlabeled.all() is ma.masked:
            break

        y_unlabeled[index_labelable] = [np.argmax(p) for p in probs[index_greater_equal]]

        X_labelable = X_unLabeled[index_unlabeled.mask]
        y_labelable = y_unlabeled[index_unlabeled.mask]

        clf.fit(X=np.append(X, X_labelable, axis=0),
                y=np.append(y, y_labelable))

    if train_is_failed:
        y_unlabeled = []
    else:
        y_unlabeled = ma.array(data=y_unlabeled, mask=index_unlabeled.mask)

    return clf, y_unlabeled

def woa_profile_from_dap(var, d, lat, lon, depth, cfg):
    """
    Monthly Climatologic Mean and Standard Deviation from WOA,
    used either for temperature or salinity.

    INPUTS
        time: [day of the year]
        lat: [-90<lat<90]
        lon: [-180<lon<180]
        depth: [meters]

    Reads the WOA Monthly Climatology NetCDF file and
    returns the corresponding WOA values of salinity or temperature mean and
    standard deviation for the given time, lat, lon, depth.
    """
    if lon < 0:
        lon = lon+360

    url = cfg['url']

    doy = int(d.strftime('%j'))
    dataset = open_url(url)

    dn = (np.abs(doy-dataset['time'][:])).argmin()
    xn = (np.abs(lon-dataset['lon'][:])).argmin()
    yn = (np.abs(lat-dataset['lat'][:])).argmin()

    if re.match("temperature\d?$", var):
        mn = ma.masked_values(dataset.t_mn.t_mn[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.t_mn.attributes['_FillValue'])
        sd = ma.masked_values(dataset.t_sd.t_sd[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.t_sd.attributes['_FillValue'])
        # se = ma.masked_values(dataset.t_se.t_se[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.t_se.attributes['_FillValue'])
        # Use this in the future. A minimum # of samples
        # dd = ma.masked_values(dataset.t_dd.t_dd[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.t_dd.attributes['_FillValue'])
    elif re.match("salinity\d?$", var):
        mn = ma.masked_values(dataset.s_mn.s_mn[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.s_mn.attributes['_FillValue'])
        sd = ma.masked_values(dataset.s_sd.s_sd[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.s_sd.attributes['_FillValue'])
        # dd = ma.masked_values(dataset.s_dd.s_dd[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.s_dd.attributes['_FillValue'])
    zwoa = ma.array(dataset.depth[:])

    ind = (depth <= zwoa.max()) & (depth >= zwoa.min())
    # Mean value profile
    f = interp1d(zwoa[~ma.getmaskarray(mn)].compressed(), mn.compressed())
    mn_interp = ma.masked_all(depth.shape)
    mn_interp[ind] = f(depth[ind])
    # The stdev profile
    f = interp1d(zwoa[~ma.getmaskarray(sd)].compressed(), sd.compressed())
    sd_interp = ma.masked_all(depth.shape)
    sd_interp[ind] = f(depth[ind])

    output = {'woa_an': mn_interp, 'woa_sd': sd_interp}

    return output

def test_no_data_available():
    """ This is a position without valid data """

    db = WOA()
    out = db['TEMP'].extract(doy=155, lat=48.1953, lon=-69.5855,
            depth=[2.0, 5.0, 6.0, 21.0, 44.0, 79.0, 5000])
    assert sorted(out.keys()) == [u't_dd', u't_mn', u't_sd', u't_se']
    for v in out:
        ma.getmaskarray(out[v]).all()

    def test_set_fields(self):
        # Tests setting fields.
        base = self.base.copy()
        mbase = base.view(mrecarray)
        mbase = mbase.copy()
        mbase.fill_value = (999999, 1e20, 'N/A')
        # Change the data, the mask should be conserved
        mbase.a._data[:] = 5
        assert_equal(mbase['a']._data, [5, 5, 5, 5, 5])
        assert_equal(mbase['a']._mask, [0, 1, 0, 0, 1])
        # Change the elements, and the mask will follow
        mbase.a = 1
        assert_equal(mbase['a']._data, [1]*5)
        assert_equal(ma.getmaskarray(mbase['a']), [0]*5)
        # Use to be _mask, now it's recordmask
        assert_equal(mbase.recordmask, [False]*5)
        assert_equal(mbase._mask.tolist(),
                     np.array([(0, 0, 0),
                               (0, 1, 1),
                               (0, 0, 0),
                               (0, 0, 0),
                               (0, 1, 1)],
                              dtype=bool))
        # Set a field to mask ........................
        mbase.c = masked
        # Use to be mask, and now it's still mask !
        assert_equal(mbase.c.mask, [1]*5)
        assert_equal(mbase.c.recordmask, [1]*5)
        assert_equal(ma.getmaskarray(mbase['c']), [1]*5)
        assert_equal(ma.getdata(mbase['c']), [asbytes('N/A')]*5)
        assert_equal(mbase._mask.tolist(),
                     np.array([(0, 0, 1),
                               (0, 1, 1),
                               (0, 0, 1),
                               (0, 0, 1),
                               (0, 1, 1)],
                              dtype=bool))
        # Set fields by slices .......................
        mbase = base.view(mrecarray).copy()
        mbase.a[3:] = 5
        assert_equal(mbase.a, [1, 2, 3, 5, 5])
        assert_equal(mbase.a._mask, [0, 1, 0, 0, 0])
        mbase.b[3:] = masked
        assert_equal(mbase.b, base['b'])
        assert_equal(mbase.b._mask, [0, 1, 0, 1, 1])
        # Set fields globally..........................
        ndtype = [('alpha', '|S1'), ('num', int)]
        data = ma.array([('a', 1), ('b', 2), ('c', 3)], dtype=ndtype)
        rdata = data.view(MaskedRecords)
        val = ma.array([10, 20, 30], mask=[1, 0, 0])

        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            rdata['num'] = val
            assert_equal(rdata.num, val)
            assert_equal(rdata.num.mask, [1, 0, 0])

 def test_set_mask(self):
     base = self.base.copy()
     mbase = base.view(mrecarray)
     # Set the mask to True .......................
     mbase.mask = masked
     assert_equal(ma.getmaskarray(mbase["b"]), [1] * 5)
     assert_equal(mbase["a"]._mask, mbase["b"]._mask)
     assert_equal(mbase["a"]._mask, mbase["c"]._mask)
     assert_equal(mbase._mask.tolist(), np.array([(1, 1, 1)] * 5, dtype=bool))
     # Delete the mask ............................
     mbase.mask = nomask
     assert_equal(ma.getmaskarray(mbase["c"]), [0] * 5)
     assert_equal(mbase._mask.tolist(), np.array([(0, 0, 0)] * 5, dtype=bool))

 def outer (self, a, b):
     "Return the function applied to the outer product of a and b."
     a = _makeMaskedArg(a)
     b = _makeMaskedArg(b)
     ma = getmask(a)
     mb = getmask(b)
     if ma is nomask and mb is nomask:
         m = None
     else:
         ma = getmaskarray(a)
         mb = getmaskarray(b)
         m = logical_or.outer(ma, mb)
     d = numpy.maximum.outer(filled(a), filled(b))
     return TransientVariable(d, mask=m)

def computeEdgeDistances(uvframe):
    """
    Create a 2D matrix @edgedists as a companion to @uvframe,
    containing for each pixel a distance to the nearest edge (more precisely,
    the nearest 0-valued pixel).

    We compute @edgedists in a floodfill fashion spreading from zero-areas
    to the middle of one-areas iteratively, with distances approximated
    on the pixel grid.

    We return a tuple (edgedists, edgedirs), where edgedirs contains information
    about the relative offset of the nearest edge piece.
    """
    # edgedists is a masked array, with only already computed values unmasked;
    # at first, uvframe == 0 already are computed (as zeros)
    edgedists = ma.array(numpy.zeros(uvframe.shape, dtype = numpy.float), mask = (uvframe > 0))
    edgedirs = ma.array(numpy.zeros(uvframe.shape, dtype = (numpy.float, 2)), mask = [[[j,j] for j in i] for i in uvframe > 0])
    #numpy.set_printoptions(threshold=numpy.nan)
    #print edgedists
    #print edgedirs

    flood_spread = scipy.ndimage.morphology.generate_binary_structure(2, 2)
    neighbor_ofs = [[-1,-1],[-1,0],[-1,1], [0,-1],[0,0],[0,1],  [1,-1],[1,0],[1,1]]
    s2 = math.sqrt(2)
    neighbor_dist = [s2,1,s2, 1,0,1, s2,1,s2]

    while ma.getmaskarray(edgedists).any():
        # scan masked area for any elements that have unmasked neighbors
        done_mask = numpy.invert(ma.getmaskarray(edgedists))
        todo_mask = done_mask ^ scipy.ndimage.binary_dilation(done_mask, flood_spread)
        #print_mask(todo_mask)
        for i in numpy.transpose(numpy.nonzero(todo_mask)):
            neighbor_val = ma.array([
                    edge_dist_if_within(edgedists, i + ofs) + dist
                        for ofs, dist in zip(neighbor_ofs, neighbor_dist)
                ])
            nearestnei = ma.argmin(neighbor_val)

            # We assert that this update never affects value other fields
            # visited later in this iteration of floodfill
            edgedists[tuple(i)] = neighbor_val[nearestnei]

            nearestneicoord = i + neighbor_ofs[nearestnei]
            #print "-", nearestneicoord, edgedirs[tuple(nearestneicoord)]
            edgedirs[tuple(i)] = edgedirs[tuple(nearestneicoord)] + tuple(neighbor_ofs[nearestnei])
            #print "+", i, edgedirs[tuple(i)]

    return (edgedists.data, edgedirs.data)

 def test_off_map(self) :
     Data = self.blocks[0]
     Data.calc_freq()
     map = self.map
     map[:,:,:] = 0.0
     Data.data[:,:,:,:] = 0.0
     # Rig the pointing but put one off the map.
     def rigged_pointing() :
         Data.ra = map.get_axis('ra')[range(10)]
         Data.dec = map.get_axis('dec')[range(10)]
         Data.ra[3] = Data.ra[3] - 8.0
     Data.calc_pointing = rigged_pointing
     smd.sub_map(Data, map)
     self.assertTrue(sp.alltrue(ma.getmaskarray(Data.data[3,:,:,:])))
     self.assertTrue(sp.alltrue(sp.logical_not(
                 ma.getmaskarray((Data.data[[0,1,2,4,5,6,7,8,9],:,:,:])))))

 def calculate_theta(self, Xm, p_y_given_x):
     """Estimate marginal parameters from data and expected latent labels."""
     theta = []
     for i in range(self.n_visible):
         not_missing = np.logical_not(ma.getmaskarray(Xm)[:, i])
         theta.append(self.estimate_parameters(Xm.data[not_missing, i], p_y_given_x[:, not_missing]))
     return np.array(theta)

    def test(self):
        self.flags = {}
        try:
            threshold = self.cfg['threshold']
        except:
            print("Deprecated cfg format. It should contain a threshold item.")
            threshold = self.cfg

        try:
            flag_good = self.cfg['flag_good']
            flag_bad = self.cfg['flag_bad']
        except:
            print("Deprecated cfg format. It should contain flag_good & flag_bad.")
            flag_good = 1
            flag_bad = 4

        assert (np.size(threshold) == 1) and \
                (threshold is not None) and \
                (np.isfinite(threshold))   

        flag = np.zeros(self.data[self.varname].shape, dtype='i1')
        flag[np.nonzero(self.features['bin_spike'] > threshold)] = flag_bad
        flag[np.nonzero(self.features['bin_spike'] <= threshold)] = flag_good
        flag[ma.getmaskarray(self.data[self.varname])] = 9
        self.flags['bin_spike'] = flag

    def start(self, observation):
        '''
        Handle the first iteration.
        
        Keyword Arguments:
        
        observation (int) The index of the most recent observation [0:num_states]
        
        Returns: 
        next_action (int) The index of the next action to take [0:num_actions]
        '''
        next_state = observation
        
        # the valid actions are at unmasked array positions, so compute the inverse
        # of the mask of the relevant row
        valid_action_flags = ~ma.getmaskarray(self._q_table[next_state,:])
        
        # get the indices of the valid actions and pick one at random
        next_action = random.choice(np.flatnonzero(valid_action_flags))

        # store off state for the next iteration
        self._last_state = next_state
        self._last_action = next_action
        
        # increment state visitation table
        self.update_visitation_table(next_state, next_action)
        self._epoch_num +=1
          
        return next_action

def wmean_bandpass_1D_serial(data, lshorterpass, llongerpass, t=None,
        method='hann', axis=0):
    """ Equivalent to wmean_1D_serial, but it is a bandpass

        Input:
            - data: np.array or ma.maked_array, nD
            - lshorterpass: The size of the highpass filter, i.e. shorter
                wavelenghts are preserved. It is in the same unit of t.
            - llongerpass: The size of the lowpass filter, i.e.longer
                wavelenghts are preserved. It is in the same unit of t.
	    - t: is the scale of the choosed axis, 1D. If not
                defined, it will be considered a sequence.
            - method: ['hann', 'hamming', 'blackman']
                Defines the weight function type
            - axis: Dimension which the filter will be applied
    """
    assert False, "There is a BUG here"

    assert axis <= data.ndim, "Invalid axis!"

    # If necessary, move the axis to be filtered for the first axis
    if axis != 0:
        data_smooth = wmean_bandpass_1D_serial(data.swapaxes(0, axis),
                lshorterpass = lshorterpass,
                llongerpass = llongerpass,
                t = t,
                method = method,
                axis = 0)

        return data_smooth.swapaxes(0, axis)
    # Below here, the filter will be always applied on axis=0

    # If t is not given, creates a regularly spaced t
    if t is None:
        print "The scale along the choosed axis weren't defined. I'll consider a constant sequence."
	t = np.arange(data.shape[axis])

    assert t.shape == (data.shape[axis],), "Invalid size of t."

    # ----
    winfunc = window_func(method)

    data_smooth = ma.masked_all(data.shape)

    if data.ndim==1:
        (I,) = np.nonzero(~ma.getmaskarray(data))
        for i in I:
            # First remove the high frequency
            tmp = _convolve_1D(t[i], t, llongerpass, winfunc, data)
            # Then remove the low frequency
            data_smooth[i] = tmp - \
                    _convolve_1D(t[i], t, lshorterpass, winfunc, tmp)

    else:
        I = data.shape[1]
        for i in range(I):
            data_smooth[:,i] = wmean_bandpass_1D_serial(data[:,i],
                    lshorterpass, llongerpass, t, method, axis)

    return data_smooth