def rectify_obs(obs):
    """Make sure the passed obs dictionary conforms to code expectations,
    and make simple fixes when possible.
    """
    k = obs.keys()
    if 'maggies' not in k:
        obs['maggies'] = None
        obs['maggies_unc'] = None
    if 'spectrum' not in k:
        obs['spectrum'] = None
        obs['unc'] = None
    if obs['maggies'] is not None:
        assert (len(obs['filters']) == len(obs['maggies']))
        assert ('maggies_unc' in k)
        assert ((len(obs['maggies']) == len(obs['maggies_unc'])) or
                (np.size(obs['maggies_unc'] == 1)))
        m = obs.get('phot_mask', np.ones(len(obs['maggies']), dtype=bool))
        obs['phot_mask'] = (m * np.isfinite(obs['maggies']) *
                            np.isfinite(obs['maggies_unc']) *
                            (obs['maggies_unc'] > 0))
        try:
            obs['filternames'] = [f.name for f in obs['filters']]
        except:
            pass

    if 'logify_spectrum' not in k:
        obs['logify_spectrum'] = False
    if obs['spectrum'] is not None:
        assert (len(obs['wavelength']) == len(obs['spectrum']))
        assert ('unc' in k)
        np.errstate(invalid='ignore')
        m = obs.get('mask', np.ones(len(obs['wavelength']), dtype=bool))
        obs['mask'] = (m * np.isfinite(obs['spectrum']) *
                       np.isfinite(obs['unc']) * (obs['unc'] > 0))
    return obs

def getEff(s, cut, comp='joint', reco=True):

    eff, sig, relerr = {},{},{}
    a = np.log10(s['MC_energy'])
    Ebins = getEbins()
    Emids = getMids(Ebins)
    erangeDict = getErange()

    c0 = cut
    if comp != 'joint':
        compcut = s['comp'] == comp
        c0 = cut * compcut

    # Set radii for finding effective area
    rDict = {}
    keys = ['low', 'mid', 'high']
    for key in keys:
        rDict[key] = np.array([600, 800, 1100, 1700, 2600, 2900])
    rDict['low'][1] = 600
    Ebreaks = np.array([4, 5, 6, 7, 8, 9])
    rgrp = np.digitize(Emids, Ebreaks) - 1

    for key in keys:

        # Get efficiency and sigma
        simcut = np.array([sim in erangeDict[key] for sim in s['sim']])
        k = np.histogram(a[c0*simcut], bins=Ebins)[0]
        #k = Nfinder(a, c0*simcut)
        n = s['MC'][comp][key].astype('float')
        eff[key], sig[key], relerr[key] = np.zeros((3, len(k)))
        with np.errstate(divide='ignore', invalid='ignore'):
            eff[key] = k / n
            var = (k+1)*(k+2)/((n+2)*(n+3)) - (k+1)**2/((n+2)**2)
        sig[key] = np.sqrt(var)

        # Multiply by throw area
        r = np.array([rDict[key][i] for i in rgrp])
        eff[key] *= np.pi*(r**2)
        sig[key] *= np.pi*(r**2)

        # Deal with parts of the arrays with no information
        for i in range(len(eff[key])):
            if n[i] == 0:
                eff[key][i] = 0
                sig[key][i] = np.inf

    # Combine low, mid, and high energy datasets
    eff_tot = (np.sum([eff[key]/sig[key] for key in keys], axis=0) /
            np.sum([1/sig[key] for key in keys], axis=0))
    sig_tot = np.sqrt(1 / np.sum([1/sig[key]**2 for key in keys], axis=0))
    with np.errstate(divide='ignore'):
        relerr  = sig_tot / eff_tot

    # UGH : find better way to do this
    if reco:
        eff_tot = eff_tot[20:]
        sig_tot = sig_tot[20:]
        relerr  = relerr[20:]

    return eff_tot, sig_tot, relerr

 def test_no_scaling(self):
     # Test writing image converting types when no scaling
     img_class = self.image_class
     hdr_class = img_class.header_class
     hdr = hdr_class()
     supported_types = supported_np_types(hdr)
     slope = 2
     inter = 10 if hdr.has_data_intercept else 0
     for in_dtype, out_dtype in itertools.product(
         FLOAT_TYPES + IUINT_TYPES,
         supported_types):
         # Need to check complex scaling
         mn_in, mx_in = _dt_min_max(in_dtype)
         arr = np.array([mn_in, -1, 0, 1, 10, mx_in], dtype=in_dtype)
         img = img_class(arr, np.eye(4), hdr)
         img.set_data_dtype(out_dtype)
         img.header.set_slope_inter(slope, inter)
         with np.errstate(invalid='ignore'):
             rt_img = bytesio_round_trip(img)
         with suppress_warnings():  # invalid mult
             back_arr = rt_img.get_data()
         exp_back = arr.copy()
         if in_dtype not in COMPLEX_TYPES:
             exp_back = arr.astype(float)
         if out_dtype in IUINT_TYPES:
             with np.errstate(invalid='ignore'):
                 exp_back = np.round(exp_back)
             exp_back = np.clip(exp_back, *shared_range(float, out_dtype))
             exp_back = exp_back.astype(out_dtype).astype(float)
         else:
             exp_back = exp_back.astype(out_dtype)
         # Allow for small differences in large numbers
         with suppress_warnings():  # invalid value
             assert_allclose_safely(back_arr,
                                    exp_back * slope + inter)

def test_3d():
    x = np.array([[9.0, 3.0, nan, nan, 9.0, nan],
                  [1.0, 1.0, 1.0, nan, nan, nan],
                  [2.0, 2.0, 0.1, nan, 1.0, nan],  # 0.0 kills geometric mean
                  [3.0, 9.0, 2.0, nan, nan, nan],
                  [4.0, 4.0, 3.0, 9.0, 2.0, nan],
                  [5.0, 5.0, 4.0, 4.0, nan, nan]])
    sectors = ['a', 'b', 'a', 'b', 'a', 'c']
    x = np.dstack((x,x))
    
    for func in funcs_one:
        xc = x.copy()
        args = (xc,)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args
        
    for func in funcs_oneint:
        xc = x.copy()
        args = (xc, 2)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args
        
    for func in funcs_onefrac:
        xc = x.copy()
        args = (xc, -1, 0.5)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args
    
    for func in funcs_sect:
        xc = x.copy()
        args = (xc, sectors)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args

    def _convertImagesToUint8(self, imageR, imageG, imageB):
        """Use the mapping to convert images imageR, imageG, and imageB to a triplet of uint8 images
        """
        imageR = imageR - self.minimum[0]  # n.b. makes copy
        imageG = imageG - self.minimum[1]
        imageB = imageB - self.minimum[2]

        fac = self.mapIntensityToUint8(self.intensity(imageR, imageG, imageB))

        imageRGB = [imageR, imageG, imageB]
        with np.errstate(invalid="ignore"):  # suppress NAN warnings
            for c in imageRGB:
                c *= fac
                # individual bands can still be < 0, even if fac isn't
                c[c < 0] = 0

        pixmax = self._uint8Max
        # copies -- could work row by row to minimise memory usage
        r0, g0, b0 = imageRGB

        # n.b. np.where can't and doesn't short-circuit
        with np.errstate(invalid='ignore', divide='ignore'):
            for i, c in enumerate(imageRGB):
                c = np.where(r0 > g0,
                             np.where(r0 > b0,
                                      np.where(r0 >= pixmax, c*pixmax/r0, c),
                                      np.where(b0 >= pixmax, c*pixmax/b0, c)),
                             np.where(g0 > b0,
                                      np.where(g0 >= pixmax, c*pixmax/g0, c),
                                      np.where(b0 >= pixmax, c*pixmax/b0, c))).astype(np.uint8)
                c[c > pixmax] = pixmax

                imageRGB[i] = c

        return imageRGB

def _solve_quadratic(a__, b__, c__, min_val=0.0, max_val=1.0):
    """Solve quadratic equation and return the valid roots from interval
    [*min_val*, *max_val*]

    """

    def int_and_float_to_numpy(val):
        if not isinstance(val, np.ndarray):
            if isinstance(val, (int, float)):
                val = [val]
            val = np.array(val)
        return val

    a__ = int_and_float_to_numpy(a__)
    b__ = int_and_float_to_numpy(b__)
    c__ = int_and_float_to_numpy(c__)

    discriminant = b__ * b__ - 4 * a__ * c__

    # Solve the quadratic polynomial
    with np.errstate(invalid='ignore', divide='ignore'):
        x_1 = (-b__ + np.sqrt(discriminant)) / (2 * a__)
        x_2 = (-b__ - np.sqrt(discriminant)) / (2 * a__)

    # Find valid solutions, ie. 0 <= t <= 1
    x__ = x_1.copy()
    with np.errstate(invalid='ignore'):
        idxs = (x_1 < min_val) | (x_1 > max_val)
    x__[idxs] = x_2[idxs]

    with np.errstate(invalid='ignore'):
        idxs = (x__ < min_val) | (x__ > max_val)
    x__[idxs] = np.nan

    return x__

 def test_scaling(self):
     # Test integer scaling from float
     # Analyze headers cannot do float-integer scaling
     hdr = self.header_class()
     assert_true(hdr.default_x_flip)
     shape = (1,2,3)
     hdr.set_data_shape(shape)
     hdr.set_data_dtype(np.float32)
     data = np.ones(shape, dtype=np.float64)
     S = BytesIO()
     # Writing to float datatype doesn't need scaling
     hdr.data_to_fileobj(data, S)
     rdata = hdr.data_from_fileobj(S)
     assert_array_almost_equal(data, rdata)
     # Now test writing to integers
     hdr.set_data_dtype(np.int32)
     # Writing to int needs scaling, and raises an error if we can't scale
     if not hdr.has_data_slope:
         assert_raises(HeaderTypeError, hdr.data_to_fileobj, data, BytesIO())
     # But if we aren't scaling, convert the floats to integers and write
     with np.errstate(invalid='ignore'):
         hdr.data_to_fileobj(data, S, rescale=False)
     rdata = hdr.data_from_fileobj(S)
     assert_true(np.allclose(data, rdata))
     # This won't work for floats that aren't close to integers
     data_p5 = data + 0.5
     with np.errstate(invalid='ignore'):
         hdr.data_to_fileobj(data_p5, S, rescale=False)
     rdata = hdr.data_from_fileobj(S)
     assert_false(np.allclose(data_p5, rdata))

def _findAndCategorizeContours(binary):
    # see auromat.utils.outline_opencv for why we need to pad the image
    imCv = np.zeros((binary.shape[0]+2, binary.shape[1]+2), dtype=np.uint8)
    imCv[1:-1,1:-1] = binary
    contours,_ = cv.findContours(imCv, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)    
    contours = np.asarray(contours)
    contours -= 1

    area = np.asarray([cv.contourArea(c) for c in contours])
    rectAxes = np.asarray([cv.minAreaRect(c)[1] for c in contours])

    # isBigContour needs to be big enough to not discard bigger stars!
    # this could leave some spacecraft structures intact which may or may not confuse astrometry    
    # TODO the ratio below should depend on the estimated celestial pixel scale
    #      and the exposure time (longer exposure = longer star trails) 
    bigContourAreaRatio = 0.000013 # 0.0013% of the image area (~160 pixels for 12MP images)
    bigContourArea = bigContourAreaRatio*(binary.shape[0]*binary.shape[1])
    isBigContour = area > int(bigContourArea)
    isSmallContour = ~isBigContour
    
    longRatioThreshold = 5
    with np.errstate(divide='ignore', invalid='ignore'): # division produces nans and infs
        rectRatio = rectAxes[:,0]/rectAxes[:,1] if len(contours) > 0 else np.array([])

    with np.errstate(invalid='ignore'):
        isLongContour = np.logical_and(area > 20, # exclude very tiny long contours (could be stars) and inf ratios
                                       np.logical_or(rectRatio > longRatioThreshold, 
                                                     rectRatio < 1/longRatioThreshold)
                                       )
    isSmallLongContour = np.logical_and(isSmallContour, isLongContour)
    isSmallShortContour = np.logical_and(isSmallContour, ~isLongContour)
    
    return contours, area, isBigContour, isSmallLongContour, isSmallShortContour

    def test_modulo(self):
        # GH3590, modulo as ints
        p = DataFrame({'first': [3, 4, 5, 8], 'second': [0, 0, 0, 3]})

        # this is technically wrong as the integer portion is coerced to float
        # ###
        expected = DataFrame({'first': Series([0, 0, 0, 0], dtype='float64'),
                              'second': Series([np.nan, np.nan, np.nan, 0])})
        result = p % p
        assert_frame_equal(result, expected)

        # numpy has a slightly different (wrong) treatement
        with np.errstate(all='ignore'):
            arr = p.values % p.values
        result2 = DataFrame(arr, index=p.index,
                            columns=p.columns, dtype='float64')
        result2.iloc[0:3, 1] = np.nan
        assert_frame_equal(result2, expected)

        result = p % 0
        expected = DataFrame(np.nan, index=p.index, columns=p.columns)
        assert_frame_equal(result, expected)

        # numpy has a slightly different (wrong) treatement
        with np.errstate(all='ignore'):
            arr = p.values.astype('float64') % 0
        result2 = DataFrame(arr, index=p.index, columns=p.columns)
        assert_frame_equal(result2, expected)

        # not commutative with series
        p = DataFrame(np.random.randn(10, 5))
        s = p[0]
        res = s % p
        res2 = p % s
        self.assertFalse(np.array_equal(res.fillna(0), res2.fillna(0)))

def test_nanvar_issue60():
    "nanvar regression test (issue #60)"

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")

        f = bn.nanvar([1.0], ddof=1)
        with np.errstate(invalid='ignore'):
            s = bn.slow.nanvar([1.0], ddof=1)
        assert_equal(f, s, err_msg="bn.nanvar([1.0], ddof=1) wrong")

        f = bn.nanvar([1], ddof=1)
        with np.errstate(invalid='ignore'):
            s = bn.slow.nanvar([1], ddof=1)
        assert_equal(f, s, err_msg="bn.nanvar([1], ddof=1) wrong")

        f = bn.nanvar([1, np.nan], ddof=1)
        with np.errstate(invalid='ignore'):
            s = bn.slow.nanvar([1, np.nan], ddof=1)
        assert_equal(f, s, err_msg="bn.nanvar([1, nan], ddof=1) wrong")

        f = bn.nanvar([[1, np.nan], [np.nan, 1]], axis=0, ddof=1)
        with np.errstate(invalid='ignore'):
            s = bn.slow.nanvar([[1, np.nan], [np.nan, 1]], axis=0, ddof=1)
        assert_equal(f, s, err_msg="issue #60 regression")

def move_unit_maker(func, arrfunc, methods):
    "Test that different mov methods give the same results on 2d input."
    arr1 = np.array([1, 2, 3, 4, 5, 6, nan, nan, 7, 8, 9])
    arr2 = np.array([[9.0, 3.0, nan, nan, 9.0, nan],
                     [1.0, 1.0, 1.0, nan, nan, nan],
                     [2.0, 2.0, 0.1, nan, 1.0, nan],
                     [3.0, 9.0, 2.0, nan, nan, nan],
                     [4.0, 4.0, 3.0, 9.0, 2.0, nan],
                     [5.0, 5.0, 4.0, 4.0, nan, nan]]) 
    arr3 = np.arange(60).reshape(3, 4, 5)
    arr4 = np.array([nan, nan, nan])
    arrs = [arr1, arr2, arr3, arr4]
    msg = '\nfunc %s | method %s | nd %d | window %d | axis %d\n'
    for arr in arrs:
        for axis in range(arr.ndim):
            for w in range(1, arr.shape[axis]):
                actual = func(arr, window=w, axis=axis, method='loop')
                for method in methods:
                    if method == 'func_loop':
                        with np.errstate(invalid='ignore'):
                            d = move_func(arrfunc, arr, window=w, axis=axis,
                                          method='loop')
                    elif method == 'func_strides':
                        with np.errstate(invalid='ignore'):
                            d = move_func(arrfunc, arr, window=w, axis=axis,
                                          method='strides')
                    else:
                        d = func(arr, window=w, axis=axis, method=method) 
                    err_msg = msg % (func.__name__, method, arr.ndim, w, axis)
                    assert_array_almost_equal(actual, d, 10, err_msg)

def test_return_array():
    "Check that functions return a numpy array or a scalar."
    
    x = np.array([[9.0, 3.0, nan, nan, 9.0, nan],
                  [1.0, 1.0, 1.0, nan, nan, nan],
                  [2.0, 2.0, 9.0, nan, 1.0, nan],
                  [3.0, 9.0, 2.0, nan, nan, nan],
                  [4.0, 4.0, 3.0, 9.0, 2.0, nan],
                  [5.0, 5.0, 4.0, 4.0, nan, nan]])
    sectors = ['a', 'b', 'a', 'b', 'a', 'c']
    
    for func in funcs_one:
        xc = x.copy()
        args = (xc,)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args
        
    for func in funcs_oneint:
        xc = x.copy()
        args = (xc, 2)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args
        
    for func in funcs_onefrac:
        xc = x.copy()
        args = (xc, -1, 0.5)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args
        
    for func in funcs_sect:
        xc = x.copy()
        args = (xc, sectors)
        with np.errstate(invalid='ignore', divide='ignore'):
            yield check_return_array, func, args

def compute_mandelbrot():
	N_max = 50
	some_threshold = 50

	grid_interval = 1000

	# construct the grid 
	x = np.linspace(-2, 1, grid_interval)
	y = np.linspace(-1.5, 1.5, grid_interval)

	c = x[:, np.newaxis] + 1j*y[np.newaxis, :]

	# do the iteration 
	z = c 
	for v in range(N_max):
		with np.errstate(all = "ignore"): # catches overflow and invalid value errors 
			z = z**2 + c 
		
	with np.errstate(all = "ignore"): # catches overflow and invalid value errors 

		# form a 2-D boolean mask 
		mask = (abs(z) < some_threshold)

	# save the result to an image
	plt.imshow(mask.T, extent = [-2, 1, -1.5, 1.5])
	plt.gray()
	plt.savefig('mandelbrot.png')

 def test_correlation_9(self):
     "farray.correlation_9"
     x = self.a1
     y = self.a2
     x2 = np.empty((2, x.shape[0], x.shape[1]))
     x2[0] = x
     x2[1] = x
     y2 = np.empty((2, y.shape[0], y.shape[1]))
     y2[0] = y
     y2[1] = y        
     with np.errstate(invalid='ignore'):
         corr = correlation(x, y, axis=-1)
     desired = np.array([nan, 1, -1, -0.5]) 
     aae(corr, desired, err_msg="aggregate of 1d tests")
     x = self.b1
     y = self.b2
     x2 = np.empty((2, x.shape[0], x.shape[1]))
     x2[0] = x
     x2[1] = x
     y2 = np.empty((2, y.shape[0], y.shape[1]))
     y2[0] = y
     y2[1] = y        
     with np.errstate(invalid='ignore'):
         corr = correlation(x, y, axis=-1)
     desired = np.array([nan, 1, -1, -0.5]) 
     aae(corr, desired, err_msg="aggregate of 1d tests")

def test_SphericalCoordinates_bounds(pyntcloud_with_rgb_and_normals):
    scalar_field = SphericalCoordinates(
        pyntcloud=pyntcloud_with_rgb_and_normals)
    scalar_field.extract_info()

    with np.errstate(divide='ignore', invalid='ignore'):
        scalar_field.compute()

    assert all(scalar_field.to_be_added["polar"] >= 0)
    assert all(scalar_field.to_be_added["polar"] <= 180)

    assert all(scalar_field.to_be_added["azimuthal"] >= -180)
    assert all(scalar_field.to_be_added["azimuthal"] <= 180)

    scalar_field = SphericalCoordinates(
        pyntcloud=pyntcloud_with_rgb_and_normals,
        degrees=False)
    scalar_field.extract_info()

    with np.errstate(divide='ignore', invalid='ignore'):
        scalar_field.compute()

    assert all(scalar_field.to_be_added["polar"] >= 0)
    assert all(scalar_field.to_be_added["polar"] <= np.pi)

    assert all(scalar_field.to_be_added["azimuthal"] >= -np.pi)
    assert all(scalar_field.to_be_added["azimuthal"] <= np.pi)

 def predict_proba(self, X):
     """
     Return the probability that the class label is +1 given X.
     :param X: An examples-by-features NumPy matrix.
     :returns: Pr[Y=1|X]
     """
     Xstd = self._standardize_inputs(X)
     Pacc = np.zeros((2, X.shape[0]))
     for y, yprob in enumerate(self._yparam):
         for i, matrix in enumerate(self._params):
             # This is a v-length vector of probabilities conditioned on y
             probabilities = matrix[y]
             # The X values are indices into this vector!  We just sum up
             # logs (instead of multiplying).  We silence divide by zero
             # errors because that's what NumPy gives when it takes the log
             # of 0.  Thankfully, np.log(0) = -inf, and np.exp(-np.inf) = 0,
             # so we're safe here!
             with np.errstate(divide='ignore'):
                 Pacc[y] += np.log(probabilities[Xstd[:, i].astype(int)])
         Pacc[y] += np.log(yprob)
     # Now, return the probability Y is positive conditioned on X!
     probs = np.exp(Pacc)
     with np.errstate(invalid='ignore'):
         # We silence invalid value warnings here, because that's what we
         # get when we do 0/0.  In this case, we haven't seen any training
         # examples with those values, and m is set to 0, so we just go
         # 50/50 and assign 0.5!
         rv = probs[1] / (probs[0] + probs[1])
     rv[np.isnan(rv)] = 0.5
     return rv

def test_a2f_min_max():
    # Check min and max thresholding of array to file
    str_io = BytesIO()
    for in_dt in (np.float32, np.int8):
        for out_dt in (np.float32, np.int8):
            arr = np.arange(4, dtype=in_dt)
            # min thresholding
            with np.errstate(invalid='ignore'):
                data_back = write_return(arr, str_io, out_dt, 0, 0, 1, 1)
            assert_array_equal(data_back, [1, 1, 2, 3])
            # max thresholding
            with np.errstate(invalid='ignore'):
                data_back = write_return(arr, str_io, out_dt, 0, 0, 1, None, 2)
            assert_array_equal(data_back, [0, 1, 2, 2])
            # min max thresholding
            data_back = write_return(arr, str_io, out_dt, 0, 0, 1, 1, 2)
            assert_array_equal(data_back, [1, 1, 2, 2])
    # Check that works OK with scaling and intercept
    arr = np.arange(4, dtype=np.float32)
    data_back = write_return(arr, str_io, np.int, 0, -1, 0.5, 1, 2)
    assert_array_equal(data_back * 0.5 - 1, [1, 1, 2, 2])
    # Even when scaling is negative
    data_back = write_return(arr, str_io, np.int, 0, 1, -0.5, 1, 2)
    assert_array_equal(data_back * -0.5 + 1, [1, 1, 2, 2])
    # Check complex numbers
    arr = np.arange(4, dtype=np.complex64) + 100j
    with suppress_warnings():  # cast to real
        data_back = write_return(arr, str_io, out_dt, 0, 0, 1, 1, 2)
    assert_array_equal(data_back, [1, 1, 2, 2])

    def test_div(self):

        # integer div, but deal with the 0's (GH 9144)
        p = DataFrame({'first': [3, 4, 5, 8], 'second': [0, 0, 0, 3]})
        result = p / p

        expected = DataFrame({'first': Series([1.0, 1.0, 1.0, 1.0]),
                              'second': Series([nan, nan, nan, 1])})
        assert_frame_equal(result, expected)

        with np.errstate(all='ignore'):
            arr = p.values.astype('float') / p.values
        result2 = DataFrame(arr, index=p.index,
                            columns=p.columns)
        assert_frame_equal(result2, expected)

        result = p / 0
        expected = DataFrame(np.inf, index=p.index, columns=p.columns)
        expected.iloc[0:3, 1] = nan
        assert_frame_equal(result, expected)

        # numpy has a slightly different (wrong) treatement
        with np.errstate(all='ignore'):
            arr = p.values.astype('float64') / 0
        result2 = DataFrame(arr, index=p.index,
                            columns=p.columns)
        assert_frame_equal(result2, expected)

        p = DataFrame(np.random.randn(10, 5))
        s = p[0]
        res = s / p
        res2 = p / s
        self.assertFalse(np.array_equal(res.fillna(0), res2.fillna(0)))

 def solve_for(self, data, inputs, chunk, pol, sideband='USB'):
     baselines = list(baseline for baseline in data.baselines if ((baseline.left in inputs) and (baseline.right in inputs)))
     corr_matrix = zeros([SWARM_CHANNELS, len(inputs), len(inputs)], dtype=complex128)
     for baseline in baselines:
         left_i = inputs.index(baseline.left)
         right_i = inputs.index(baseline.right)
         baseline_data = data[baseline, sideband]
         complex_data = baseline_data[0::2] + 1j * baseline_data[1::2]
         corr_matrix[:, left_i, right_i] = complex_data
         corr_matrix[:, right_i, left_i] = complex_data.conj()
     this_reference = SwarmInput(self.reference.ant, chunk, pol)
     referenced_solver = partial(solve_cgains, ref=inputs.index(this_reference))
     if self.normed:
         with errstate(invalid='ignore'):
             corr_matrix = complex_nan_to_num(corr_matrix / abs(corr_matrix))
     if not self.single_chan:
         full_spec_gains = array(self.map(referenced_solver, corr_matrix))
         delays, phases = solve_delay_phase(full_spec_gains)
     else:
         full_spec_gains = array(self.map(referenced_solver, [corr_matrix.mean(axis=0),]))
         phases = (180.0/pi) * angle(full_spec_gains.mean(axis=0))
         delays = zeros(len(inputs))
     amplitudes = abs(full_spec_gains).mean(axis=0)
     with errstate(invalid='ignore'):
         efficiency = (abs(full_spec_gains.sum(axis=1)) / abs(full_spec_gains).sum(axis=1)).real
     return efficiency, vstack([amplitudes, delays, phases])

def test_nans():
    apparent_zenith = np.array([10, np.nan, 10])
    apparent_azimuth = np.array([180, 180, np.nan])
    with np.errstate(invalid='ignore'):
        tracker_data = tracking.singleaxis(apparent_zenith, apparent_azimuth,
                                           axis_tilt=0, axis_azimuth=0,
                                           max_angle=90, backtrack=True,
                                           gcr=2.0/7.0)
    expect = {'tracker_theta': np.array([0, nan, nan]),
              'aoi': np.array([10, nan, nan]),
              'surface_azimuth': np.array([90, nan, nan]),
              'surface_tilt': np.array([0, nan, nan])}
    for k, v in expect.items():
        assert_allclose(tracker_data[k], v)

    # repeat with Series because nans can differ
    apparent_zenith = pd.Series(apparent_zenith)
    apparent_azimuth = pd.Series(apparent_azimuth)
    with np.errstate(invalid='ignore'):
        tracker_data = tracking.singleaxis(apparent_zenith, apparent_azimuth,
                                           axis_tilt=0, axis_azimuth=0,
                                           max_angle=90, backtrack=True,
                                           gcr=2.0/7.0)
    expect = pd.DataFrame(np.array(
        [[ 0., 10., 90.,  0.],
         [nan, nan, nan, nan],
         [nan, nan, nan, nan]]),
        columns=['tracker_theta', 'aoi', 'surface_azimuth', 'surface_tilt'])
    assert_frame_equal(tracker_data, expect)