def assertDataAlmostEqual(self, data, reference_filename, **kwargs):
     reference_path = self.get_result_path(reference_filename)
     if self._check_reference_file(reference_path):
         kwargs.setdefault('err_msg', 'Reference file %s' % reference_path)
         with open(reference_path, 'r') as reference_file:
             stats = json.load(reference_file)
             self.assertEqual(stats.get('shape', []), list(data.shape))
             self.assertEqual(stats.get('masked', False),
                              ma.is_masked(data))
             nstats = np.array((stats.get('mean', 0.), stats.get('std', 0.),
                                stats.get('max', 0.), stats.get('min', 0.)),
                               dtype=np.float_)
             if math.isnan(stats.get('mean', 0.)):
                 self.assertTrue(math.isnan(data.mean()))
             else:
                 data_stats = np.array((data.mean(), data.std(),
                                        data.max(), data.min()),
                                       dtype=np.float_)
                 self.assertArrayAllClose(nstats, data_stats, **kwargs)
     else:
         self._ensure_folder(reference_path)
         stats = collections.OrderedDict([
             ('std', np.float_(data.std())),
             ('min', np.float_(data.min())),
             ('max', np.float_(data.max())),
             ('shape', data.shape),
             ('masked', ma.is_masked(data)),
             ('mean', np.float_(data.mean()))])
         with open(reference_path, 'w') as reference_file:
             reference_file.write(json.dumps(stats))

def sam2mat_main(args):
    region_pattern = r'^[^:]+(?::\d+-\d+)?(?:,[^:]+(?::\d+-\d+)?)?$'
    if args.region is not None and re.search(region_pattern, args.region):
        regions = args.region
    elif args.reglist is not None:
        with open(args.reglist) as f:
            regions = [line.rstrip() for line in f]
    else:
        regions = None

    if args.insam is None:
        sam_fh = sys.stdin
    else:
        sam_fh = open(args.insam, 'r')

    bdata = BinnedData(args.fai, regions=regions, resolution=args.resolution)
    bdata.read_sam(sam_fh)
    sam_fh.close()

    if args.clean:
        bdata.clean()
    if args.ice:
        bdata.iterative_correction()

    margins = bdata.dat.sum(axis=0)
    #print(margins)
    #sys.exit()

    try:
        os.makedirs(args.outdir)
    except OSError as e:
        if e.errno != errno.EEXIST:
            raise(e)
    bin_outfile = os.path.join(args.outdir, 'bins.txt.gz')
    contact_outfile = os.path.join(args.outdir, 'contacts.txt.gz')
    matrix_outfile = os.path.join(args.outdir, 'matrix.txt.gz')
    bin_f = gzip.open(bin_outfile, 'wb')
    contact_f = gzip.open(contact_outfile, 'wb')
    matrix_f = gzip.open(matrix_outfile, 'wb')

    for i,chrom1,b1 in bdata.iter_bins():
        bin_mid1 = (b1[0]+b1[1])/2
        if ma.is_masked(margins[i]):
            margin = 0
        else:
            margin = int(margins[i])
        print('{}\t{}\t{}\t{}\t{}'.format(chrom1,0,bin_mid1,margin,int(margin>0)), file=bin_f)
        if bdata.cleaned:
            print('\t'.join(bdata.dat.data[i].astype(str)), file=matrix_f)
        else:
            print('\t'.join(bdata.dat[i].astype(str)), file=matrix_f)
        for j,chrom2,b2 in bdata.iter_bins():
            bin_mid2 = (b2[0]+b2[1])/2
            contact = bdata.dat[i,j]
            if j>i and not ma.is_masked(contact) and contact > 0:
                print('{}\t{}\t{}\t{}\t{}'.format(chrom1,bin_mid1,chrom2,bin_mid2,int(contact)), file=contact_f)

    bin_f.close()
    contact_f.close()
    matrix_f.close()

def test_ephemerides_query(patch_request):
    # check values of Ceres for a given epoch
    # orbital uncertainty of Ceres is basically zero
    res = jplhorizons.Horizons(id='Ceres', location='500',
                               epochs=2451545.5).ephemerides()[0]

    assert res['targetname'] == "1 Ceres"
    assert res['datetime_str'] == "2000-Jan-01 00:00:00.000"
    assert res['solar_presence'] == ""
    assert res['flags'] == ""
    assert res['elongFlag'] == '/L'

    assert is_masked(res['AZ'])
    assert is_masked(res['EL'])
    assert is_masked(res['airmass'])
    assert is_masked(res['magextinct'])

    npt.assert_allclose(
        [2451544.5,
         188.70280, 9.09829, 34.40955, -2.68358,
         8.27, 6.83, 96.171,
         161.3828, 10.4528, 2.551099014238, 0.1744491,
         2.26315116146176, -21.9390511, 18.822054,
         95.3996, 22.5698, 292.551, 296.850,
         184.3426220, 11.7996521, 289.864329, 71.545655,
         0, 0],
        [res['datetime_jd'],
         res['RA'], res['DEC'], res['RA_rate'], res['DEC_rate'],
         res['V'], res['surfbright'], res['illumination'],
         res['EclLon'], res['EclLat'], res['r'], res['r_rate'],
         res['delta'], res['delta_rate'], res['lighttime'],
         res['elong'], res['alpha'], res['sunTargetPA'], res['velocityPA'],
         res['ObsEclLon'], res['ObsEclLat'], res['GlxLon'], res['GlxLat'],
         res['RA_3sigma'], res['DEC_3sigma']])

 def test_addTraceWithGap(self):
     """
     Tests __add__ method of the Trace class.
     """
     # set up
     tr1 = Trace(data=np.arange(1000))
     tr1.stats.sampling_rate = 200
     start = UTCDateTime(2000, 1, 1, 0, 0, 0, 0)
     tr1.stats.starttime = start
     tr2 = Trace(data=np.arange(0, 1000)[::-1])
     tr2.stats.sampling_rate = 200
     tr2.stats.starttime = start + 10
     # verify
     tr1.verify()
     tr2.verify()
     # add
     trace = tr1 + tr2
     # stats
     self.assertEquals(trace.stats.starttime, start)
     self.assertEquals(trace.stats.endtime, start + 14.995)
     self.assertEquals(trace.stats.sampling_rate, 200)
     self.assertEquals(trace.stats.npts, 3000)
     # data
     self.assertEquals(len(trace), 3000)
     self.assertEquals(trace[0], 0)
     self.assertEquals(trace[999], 999)
     self.assertTrue(is_masked(trace[1000]))
     self.assertTrue(is_masked(trace[1999]))
     self.assertEquals(trace[2000], 999)
     self.assertEquals(trace[2999], 0)
     # verify
     trace.verify()

def ham6_nearest(pixel, palette, last_color=None):
    if pixel is None or ma.is_masked(pixel):
        return ma.masked, ma.masked

    min_dist = None
    best_index = ma.masked
    best_color = ma.masked

    for i, c in enumerate(palette[:16]):
        d = color_distance(pixel, c)
        if min_dist is None or d < min_dist:
            if d == 0:
                return i, c
            min_dist = d
            best_index = i
            best_color = c

    if last_color is None or ma.is_masked(last_color):
        return best_index, best_color

    c = last_color.copy()

    for i in range(16):
        c[2] = i * 0x11
        d = color_distance(pixel, c)
        if d < min_dist:
            if d == 0:
                return i + 0x10, c
            min_dist = d
            best_index = i + 0x10
            best_color = c.copy()

    c = last_color.copy()

    for i in range(16):
        c[0] = i * 0x11
        d = color_distance(pixel, c)
        if d < min_dist:
            if d == 0:
                return i + 0x20, c
            min_dist = d
            best_index = i + 0x20
            best_color = c.copy()

    c = last_color.copy()

    for i in range(16):
        c[1] = i * 0x11
        d = color_distance(pixel, c)
        if d < min_dist:
            if d == 0:
                return i + 0x30, c
            min_dist = d
            best_index = i + 0x30
            best_color = c.copy()

    return best_index, best_color

    def add_average_field(self,
                          field_to_avg,
                          average_func=np.ma.mean,
                          valid_range=None #[20, 300]
                          ):
        """
        Will run an function over the elements of a field to reduce them
        to a single metric for each element, and add this reduced data (e.g.
        the mean value) as a new field to the DataMat.
        
        Obvious example is to compute the average pupil size in a single trial.
        
        Will honor span start and end indices if they are not masked. 
        
        Parameters:
         field_to_avg : string
             the name of the field to process
         average_func : function pointer
             a pointer to the function to use for each element
         valid_range : 2-element sequence (tuple or list)
             if not None, then minimum and maximum dictating the
             range, outside of which the data will be ignored. The data outside
             this range will be masked prior to the averaging.
        
        """
        for fieldname in [field_to_avg]:
            if fieldname not in self.fieldnames():
                raise ValueError("Required field '%s' not in Datamat." % (
                            fieldname))
        avg = []
        for dmi in self:
            dat = dmi.field(field_to_avg)[0]
            if dat is not None:
                sidx = dmi.span_start_idx[0]
                if ma.is_masked(sidx):
                    sidx = 0
                eidx = dmi.span_end_idx[0]
                if ma.is_masked(eidx):
                    eidx = -1
                spandat = dat[sidx:eidx]
                if valid_range is not None:
                    valdat = spandat[(spandat > valid_range[0]) & (spandat < valid_range[1])]
                else:
                    valdat = spandat
                datavg = average_func(valdat) if len(valdat) > 0 else np.NaN
                avg.append(datavg)
            else:
                avg.append(np.NaN)

        avg = ma.masked_invalid(avg)
        avg.fill_value = np.NaN
        fname = get_short_function_name(average_func)
        new_field = (fname) + "_" + field_to_avg

        self.add_field(new_field, avg)

def _next_non_masked_element(a, idx):
    """Return the next non masked element of a masked array.

    If an array is masked, return the next non-masked element (if the given index is masked).
    If no other unmasked points are after the given masked point, returns none.

    Parameters
    ----------
    a : array-like
        1-dimensional array of numeric values
    idx : integer
        index of requested element

    Returns
    -------
        Index of next non-masked element and next non-masked element

    """
    try:
        next_idx = idx + a[idx:].mask.argmin()
        if ma.is_masked(a[next_idx]):
            return None, None
        else:
            return next_idx, a[next_idx]
    except (AttributeError, TypeError, IndexError):
        return idx, a[idx]

def _math_op_common(cube, operation_function, new_unit, new_dtype=None,
                    in_place=False):
    _assert_is_cube(cube)

    if in_place:
        new_cube = cube
        if cube.has_lazy_data():
            new_cube.data = operation_function(cube.lazy_data())
        else:
            try:
                operation_function(cube.data, out=cube.data)
            except TypeError:
                # Non ufunc function
                operation_function(cube.data)
    else:
        new_cube = cube.copy(data=operation_function(cube.core_data()))

    # If the result of the operation is scalar and masked, we need to fix up
    # the dtype
    if new_dtype is not None \
            and not new_cube.has_lazy_data() \
            and new_cube.data.shape == () \
            and ma.is_masked(new_cube.data):
        new_cube.data = ma.masked_array(0, 1, dtype=new_dtype)

    iris.analysis.clear_phenomenon_identity(new_cube)
    new_cube.units = new_unit
    return new_cube

    def make_gene_map_2(self):
        """
        The method that takes the attributes from the array and uses
        them to create a gene map for the array.
        The gene map is a dictionary which has a binary string as a key.
        The binary string is created by creating a binary bit string of 
        an appropriate length.
        The length is calculated 
        """
        count = 0
        self.iterator_one = itertools.product(range(self.time_len), range(self.lat_len), range(self.lon_len))
        ### Assign a binary string a location and a value from the data
        print "\n"
        print "Creating gene-map dictionary... \n"
        print "Assigning valid locations to binary strings! \n"
        for x_valid in self.iterator_one:
            binary_string = bin(count)[2:]
            while len(binary_string) < self.string_length:  # removed minus one (-1) NB
                binary_string = "0" + binary_string
            self.gene_map[binary_string] = {}
            if ma.is_masked(self.array[x_valid]):
                pass
            else:
                self.gene_map[binary_string]["coordinate"] = tuple(x_valid)
                self.gene_map[binary_string]["value"] = self.array[x_valid]
                self.location_dict[x_valid[1:3]] = []
                self.location_dict_stdevs[x_valid[1:3]] = 0
                count += 1
        self.last_valid_binary_string = binary_string
        binary_string_old = binary_string
        not_valid_first = int(binary_string, 2) + 1
        not_valid_last = int("1" * (self.string_length), 2)  # added minus one just for nonmasked version NB
        self.count = count

        if self.count == self.count_non_masked:
            print "The counter corresponds with the non-masked count! \n"
        ### Pad the dictionary to give binary strings some value
        print "Assigning left over binary strings to non-existant locations! \n"
        self.iterator_two = itertools.product(range(self.time_len), range(self.lat_len), range(self.lon_len))
        # ~ for x_not_valid in range(not_valid_first, not_valid_last+1):
        count_2 = not_valid_first
        for x_not_valid in self.iterator_two:
            # ~ binary_string = bin(x_not_valid)[2:] # DOES IT NEED TO BE PADDED
            binary_string = bin(count_2)[2:]
            while len(binary_string) < self.string_length:  # removed minus one (-1) NB
                binary_string = "0" + binary_string
            self.gene_map[binary_string] = {}
            self.gene_map[binary_string]["coordinate"] = (999, 999, 999)  # x_not_valid
            self.gene_map[binary_string]["value"] = 1e09  # self.array[x_valid]
            if count_2 == not_valid_last:
                break
            else:
                count_2 += 1
        print "There are %d valid locations. \n" % count
        print "The last binary string is: ", binary_string
        print "The last binary string assigned to a valid locations is :", binary_string_old
        print "The length of binary string is: ", self.string_length
        print "The non-valid locations fall between %d and %d. \n" % (not_valid_first, not_valid_last)
        print "Is the array masked?: \n", ma.isMA(self.array)
        print "The gene-map has been created! \n"

def from_ham6(ham6, palette, background=None):
    if background is None:
        background = ma.masked
    elif isinstance(background, numbers.Integral):
        background = palette[background]

    if ma.is_masked(background) or ma.isMaskedArray(ham6):
        rgb8 = ma.empty(ham6.shape[:2] + (3,), dtype=np.uint8)
    else:
        rgb8 = np.empty(ham6.shape[:2] + (3,), dtype=np.uint8)

    for y in range(rgb8.shape[0]):
        c = background
        for x in range(rgb8.shape[1]):
            i = ham6[y, x]
            if i is ma.masked:
                ham6[y, x] = ma.masked
                continue

            if i < 0x10:
                c = palette[i]
            else:
                c = c.copy()
                c[(None, 2, 0, 1)[i >> 4]] = (i & 0xF) * 0x11

            rgb8[y, x] = c

    return rgb8

def to_ham6(img, palette, background=None, out=None):
    _debug_array(img)

    if background is None:
        background = ma.masked
    elif isinstance(background, numbers.Integral):
        background = palette[background]

    if not ma.is_masked(background) and ma.isMaskedArray(img):
        img = img.filled(background)

    if ma.isMaskedArray(img):
        ham6 = ma.empty(img.shape[:2], dtype=np.uint8)
    else:
        ham6 = np.empty(img.shape[:2], dtype=np.uint8)

    for y in range(img.shape[0]):
        c = background
        for x in range(img.shape[1]):
            i, c = ham6_nearest(img[y, x], palette, c)
            ham6[y, x] = i
            if out is not None:
                out[y, x] = c

    _debug_array(ham6)
    return ham6

def array_masked_to_nans(array):
    """
    Convert a masked array to a NumPy `ndarray` filled with NaN values. Input
    NumPy arrays with no mask are returned unchanged.
    This is used for dask integration, as dask does not support masked arrays.

    Args:

    * array:
        A NumPy `ndarray` or masked array.

    Returns:
        A NumPy `ndarray`. This is the input array if unmasked, or an array
        of floating-point values with NaN values where the mask was `True` if
        the input array is masked.

    .. note::
        The fill value and mask of the input masked array will be lost.

    .. note::
        Integer masked arrays are cast to 8-byte floats because NaN is a
        floating-point value.

    """
    if not ma.isMaskedArray(array):
        result = array
    else:
        if ma.is_masked(array):
            mask = array.mask
            new_dtype = nan_array_type(array.data.dtype)
            result = array.data.astype(new_dtype)
            result[mask] = np.nan
        else:
            result = array.data
    return result

 def err(coef):
     if ma.is_masked(v):
         res = v.flatten() - EllipticGaussian(x,y,coef).flatten()
         return res.data[res.mask==False]
     else:
         res = v.flatten() - EllipticGaussian(x,y,coef).flatten()
         return res

 def createMinMaxList(self):
     """
     Creates the minmax list. The method used is not fully accurate but the results
     should be ok.
     """
     pixel = self.win.detail
     data = self.stream.slice(self.starttime, self.endtime)[0].data 
     # Reshape and calculate point to point differences.
     per_pixel = int(len(data)//pixel)
     ptp = data[:pixel * per_pixel].reshape(pixel, per_pixel).ptp(axis=1)
     # Last pixel.
     last_pixel = data[pixel * per_pixel:]
     if len(last_pixel):
         last_pixel = last_pixel.ptp()
         if ptp[-1] < last_pixel:
             ptp[-1] = last_pixel
     self.ptp = ptp.astype('float32')
     # Create a logarithmic axis.
     if self.win.log_scale:
         self.ptp += 1
         self.ptp = np.log(self.ptp)/np.log(self.win.log_scale)
     # Make it go from 0 to 100.
     self.ptp *= 100.0/self.ptp.max()
     # Set masked arrays to zero.
     if is_masked(self.ptp):
         self.ptp.fill_value = 0.0
         self.ptp = self.ptp.filled()
     # Assure that very small values are also visible. Only true gaps are 0
     # and will stay 0.
     self.ptp[(self.ptp > 0) & (self.ptp < 0.5)] = 0.5

def load_ham6(fname, palette, background=None, mtimes=None, format=None):
    stem, __ = os.path.splitext(fname)
    fname_cache = stem + '.cache'

    masked = background is None or ma.is_masked(background)

    if mtimes is None:
        mtimes = ()
    elif not isinstance(mtimes, collections.Iterable):
        mtimes = (mtimes,)

    try:
        if os.path.exists(fname):
            mtimes = itertools.chain(mtimes, (os.path.getmtime(fname),))

        for mtime in mtimes:
            if os.path.getmtime(fname_cache) < mtime:
                raise OSError("cache file out of date")

        ham6 = np.genfromtxt(fname_cache, dtype=np.uint8,
                             missing_values=masked and '--',
                             usemask=masked, loose=False, invalid_raise=True)
    except OSError:
        rgb = load_image(fname, masked=masked, format=format)
        ham6 = to_ham6(rgb, palette, background=background)
        out = np.array(ham6, dtype=np.str_)
        if masked:
            out[ma.getmaskarray(ham6)] = '--'
        np.savetxt(fname_cache, out, fmt='%2s', delimiter=' ')

    return ham6

 def __call__(self,array):
     masked = ma.is_masked(array)
     if self.method is 'basemap':
        return basemap.interp(array, self.xin, self.yin, self.xout, self.yout, checkbounds=False, masked=masked, order=1)
     elif self.method is 'scipy':
        import scipy.interpolate
        interp = scipy.interpolate.interp2d(self.xin, self.yin, array, kind='linear')
        a1d = interp(self.xout[0,:],self.yout[:,0])
        return npy.reshape(a1d,self.yout.shape)

 def feature_from_tile(self,tile,out):
     if ma.is_masked(tile):
         tile = tile.compressed()
         
     tile = tile.reshape(-1)
     
     out[:] = np.bincount(tile,minlength=self.bin_count)
     
     out/=np.sum(out)

def unwrap_py(inph,in_p=(), uv=2*pi):
    """Return the input matrix unwraped the valu given in uv
    
    The same as unwrapv, but using for-s, written in python
    """
    if not is_masked(inph):
        fasei=MaskedArray(inph, isnan(inph))
    else:
        fasei=inph
        
    nx, ny=(fasei.shape[0],fasei.shape[1]) 
    
    # If the initial unwraping point is not given, take the center of the image
    # as initial coordinate
    if in_p==():
        in_p=(int(nx/2),int(ny/2))

    # Create a temporal space to mark if the points are already unwrapped
    # 0 the point has not been unwrapped
    # 1 the point has not been unwrapped, but it is in the unwrapping list 
    # 2 the point was already unwrapped

    fl=zeros((nx, ny))

    # List containing the points to unwrap
    l_un=[in_p]
    fl[in_p]=1

    # unwrapped values
    faseo=fasei.copy()    
    
    while len(l_un)>0:
        # remove the first value from the list
        cx, cy=l_un.pop(0)
    
        # Put the coordinates of unwrapped the neigbors in the list
        # And check for wrapping
        nv=0
        wv=0    
        
        for i in range(cx-1, cx+2):
            for j in range(cy-1, cy+2):
                if (i>-1) and (i<nx) and (j>-1) and (j<ny):
                    if (fl[i, j]==0)&(faseo.mask[i, j]==False):
                        fl[i, j]=1
                        l_un.append((i, j))
                    elif fl[i, j]==2:
                        wv=wv+rint((faseo[i, j]-faseo[cx, cy])/uv)
                        nv=nv+1        
        if nv!=0: 
            wv=wv/nv

        fl[cx, cy]=2
        faseo[cx, cy]=faseo[cx, cy]+wv*uv
        
    return faseo

    def _numpy_interpolation(self, point_num, eval_points):
        """

        Parameters
        ----------
        point_num: int
            Index of class position in values list
        eval_points: ndarray
            Inputs used to evaluate class member function

        Returns
        -------
        ndarray: output from member function
        """
        is_masked = ma.is_masked(eval_points)

        shape = point_num.shape
        ev_shape = eval_points.shape

        vals = self.values[point_num.ravel()]
        eval_points = np.repeat(eval_points, shape[1], axis=0)
        it = np.arange(eval_points.shape[0])

        it = np.repeat(it, eval_points.shape[1], axis=0)

        eval_points = eval_points.reshape(
            eval_points.shape[0] * eval_points.shape[1],
            eval_points.shape[-1]
        )

        scaled_points = eval_points.T
        if is_masked:
            mask = np.invert(ma.getmask(scaled_points[0]))
        else:
            mask = np.ones_like(scaled_points[0], dtype=bool)

        it = ma.masked_array(it, mask)
        scaled_points[0] = (
            (scaled_points[0] - (self._bounds[0][0])) /
            (self._bounds[0][1] - self._bounds[0][0])
        ) * (vals.shape[-2] - 1)
        scaled_points[1] += (
            (scaled_points[1] - (self._bounds[1][0])) /
            (self._bounds[1][1] - self._bounds[1][0])
        ) * (vals.shape[-1] - 1)
        scaled_points = np.vstack((it, scaled_points))

        output = np.zeros(scaled_points.T.shape[:-1])
        output[mask] = map_coordinates(vals, scaled_points.T[mask].T, order=1)

        new_shape = (*shape, ev_shape[-2])
        output = output.reshape(new_shape)

        return ma.masked_array(output, mask=mask)

def ssh2psi(lon,lat,ssh):
    from numpy import pi,sin,ndim,ma
    from pylab import meshgrid
    g = 9.8; #r = 6371.e3
    omega = 0.729e-4
    if ndim(lon)==1:
        lon,lat = meshgrid(lon,lat)
    f=2.0*omega*sin(lat*pi/180.0)
    psi = g/f * ssh
    if ma.is_masked(ssh):
        psi=ma.array(psi,mask=ssh.mask)
    return psi