def plot(self):
        ax = self.subplot

        xs = planets["Pl. Sol"]
        ys = planets["Pl. Grav"]
        self.active = w = logical_and(isfinite(xs), isfinite(ys)).nonzero()[0]
        xs, ys, clrs = xs[w], ys[w], colors[w]

        ax.set_xscale("log")
        ax.set_yscale("log")

        if not self.issmall:
            artist = ax.scatter(xs, ys, s=80, c=clrs, lw=1, picker=80.0, zorder=1000)
            ax.grid(True, which="both", ls="-", c="#222222")
            ax.axvspan(1.0 / 4, 1.0 / 0.56, facecolor="#111111", zorder=-1000)
            ax.set_xlabel(_(r"Sunlight strength (Earth units)"))
            ax.set_ylabel(_(r"Gravity strength (Earth units)"))
            # ax.xaxis.set_major_formatter(FormatStrFormatter(r'$%0.2f$'))
            # ax.yaxis.set_major_formatter(FormatStrFormatter(r'$%0.2f$'))
            ax.xaxis.set_major_formatter(EPFormatter(True))
            ax.yaxis.set_major_formatter(EPFormatter())
        else:
            artist = ax.scatter(xs, ys, s=40, c=clrs, lw=1, picker=80.0)
            ax.set_xlabel(_(r"Sunlight vs Gravity"))
            ax.xaxis.set_major_locator(NullLocator())
            ax.yaxis.set_major_locator(NullLocator())

        # ax.axis('scaled')
        ax.set_xlim(1e-9, 2e4)
        ax.set_ylim(1e-1, 1e2)

        return artist

def measure_image_moments(image):
    """Compute 0th, 1st and 2nd moments of an image.

    NaN values are ignored in the computation.

    Parameters
    ----------
    image :`astropy.io.fits.ImageHDU`
        Image to measure on.

    Returns
    -------
    image moments : list
        List of image moments:
        [A, x_cms, y_cms, x_sigma, y_sigma, sqrt(x_sigma * y_sigma)]
    """
    x, y = coordinates(image, lon_sym=True)
    A = image.data[np.isfinite(image.data)].sum()

    # Center of mass
    x_cms = (x * image.data)[np.isfinite(image.data)].sum() / A
    y_cms = (y * image.data)[np.isfinite(image.data)].sum() / A

    # Second moments
    x_var = ((x - x_cms) ** 2 * image.data)[np.isfinite(image.data)].sum() / A
    y_var = ((y - y_cms) ** 2 * image.data)[np.isfinite(image.data)].sum() / A
    x_sigma = np.sqrt(x_var)
    y_sigma = np.sqrt(y_var)

    return A, x_cms, y_cms, x_sigma, y_sigma, np.sqrt(x_sigma * y_sigma)

def censor_diagnosis(genotype_file,phenotype_file,final_pfile, final_gfile,field ='na',start_time=float('nan'),end_time=float('nan')):
        import pandas as pd
        import numpy as np
        genotypes = pd.read_csv(genotype_file)
        phenotypes = pd.read_csv(phenotype_file)
        mg=pd.merge(phenotypes,genotypes,on='id')
        if np.isnan(start_time) and np.isnan(end_time):
                print("Choose appropriate time period")
        if field=='na':
                if np.isfinite(start_time) and np.isnan(end_time):
                        final = mg[mg['AgeAtICD']>=start_time]
                elif np.isnan(start_time) and np.isfinite(end_time):
                        final = mg[mg['AgeAtICD']<=end_time]
                else:
                        final = mg[(mg['AgeAtICD']>=start_time)&(mg['AgeAtICD']<=end_time)]

        else:
                mg['diff']=mg[field]-mg['AgeAtICD']
                if np.isfinite(start_time) and np.isnan(end_time):
                        final = mg[(mg['diff']>=start_time)|(np.isnan(mg['diff']))]
                elif np.isnan(start_time) and np.isfinite(end_time):
                        final = mg[(mg['diff']<=end_time)|(np.isnan(mg['diff']))]
                else:
                        final = mg[(mg['diff']>=start_time)&(mg['diff']<=end_time)|(np.isnan(mg['diff']))]
        final[['id','icd9','AgeAtICD']].to_csv(final_pfile)
        final_gp = final.drop_duplicates('id')
        del final_gp['icd9']
        del final_gp['AgeAtICD']
        final_gp.to_csv(final_gfile)

def _assert_all_finite(X):
    """Like assert_all_finite, but only for ndarray."""
    X = np.asanyarray(X)
    if (X.dtype.char in np.typecodes['AllFloat'] and not np.isfinite(X.sum())
            and not np.isfinite(X).all()):
        raise ValueError("Input contains NaN, infinity"
                         " or a value too large for %r." % X.dtype)

def _plot(x, mph, mpd, threshold, edge, valley, ax, ind):
    """Plot results of the detect_peaks function, see its help."""
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print('matplotlib is not available.')
    else:
        if ax is None:
            _, ax = plt.subplots(1, 1, figsize=(8, 4))

        ax.plot(x, 'b', lw=1)
        if ind.size:
            label = 'valley' if valley else 'peak'
            label = label + 's' if ind.size > 1 else label
            ax.plot(ind, x[ind], '+', mfc=None, mec='r', mew=2, ms=8,
                    label='%d %s' % (ind.size, label))
            ax.legend(loc='best', framealpha=.5, numpoints=1)
        ax.set_xlim(-.02*x.size, x.size*1.02-1)
        ymin, ymax = x[np.isfinite(x)].min(), x[np.isfinite(x)].max()
        yrange = ymax - ymin if ymax > ymin else 1
        ax.set_ylim(ymin - 0.1*yrange, ymax + 0.1*yrange)
        ax.set_xlabel('Data #', fontsize=14)
        ax.set_ylabel('Amplitude', fontsize=14)
        mode = 'Valley detection' if valley else 'Peak detection'
        ax.set_title("%s (mph=%s, mpd=%d, threshold=%s, edge='%s')"
                     % (mode, str(mph), mpd, str(threshold), edge))
        # plt.grid()
        plt.show()

    def isclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False):
        def within_tol(x, y, atol, rtol):
            result = np.less_equal(np.abs(x-y), atol + rtol * np.abs(y))
            if np.isscalar(a) and np.isscalar(b):
                result = np.bool(result)
            return result

        x = np.array(a, copy=False, subok=True, ndmin=1)
        y = np.array(b, copy=False, subok=True, ndmin=1)
        xfin = np.isfinite(x)
        yfin = np.isfinite(y)
        if np.all(xfin) and np.all(yfin):
            return within_tol(x, y, atol, rtol)
        else:
            finite = xfin & yfin
            cond = np.zeros_like(finite, subok=True)
            # Because we're using boolean indexing, x & y must be the same shape.
            # Ideally, we'd just do x, y = broadcast_arrays(x, y). It's in
            # lib.stride_tricks, though, so we can't import it here.
            x = x * np.ones_like(cond)
            y = y * np.ones_like(cond)
            # Avoid subtraction with infinite/nan values...
            cond[finite] = within_tol(x[finite], y[finite], atol, rtol)
            # Check for equality of infinite values...
            cond[~finite] = (x[~finite] == y[~finite])
            if equal_nan:
                # Make NaN == NaN
                cond[np.isnan(x) & np.isnan(y)] = True
            return cond

    def _set_minmax(self):
        data = self._get_fast_data()
        try:
            self.maxval = numpy.nanmax(data)
            self.minval = numpy.nanmin(data)
        except Exception:
            self.maxval = 0
            self.minval = 0

        # TODO: see if there is a faster way to ignore infinity
        try:
            if numpy.isfinite(self.maxval):
                self.maxval_noinf = self.maxval
            else:
                self.maxval_noinf = numpy.nanmax(data[numpy.isfinite(data)])
        except:
            self.maxval_noinf = self.maxval

        try:
            if numpy.isfinite(self.minval):
                self.minval_noinf = self.minval
            else:
                self.minval_noinf = numpy.nanmin(data[numpy.isfinite(data)])
        except:
            self.minval_noinf = self.minval

def test_underflow_or_overlow():
    with np.errstate(all="raise"):
        # Generate some weird data with hugely unscaled features
        rng = np.random.RandomState(0)
        n_samples = 100
        n_features = 10

        X = rng.normal(size=(n_samples, n_features))
        X[:, :2] *= 1e300
        assert_true(np.isfinite(X).all())

        # Use MinMaxScaler to scale the data without introducing a numerical
        # instability (computing the standard deviation naively is not possible
        # on this data)
        X_scaled = MinMaxScaler().fit_transform(X)
        assert_true(np.isfinite(X_scaled).all())

        # Define a ground truth on the scaled data
        ground_truth = rng.normal(size=n_features)
        y = (np.dot(X_scaled, ground_truth) > 0.0).astype(np.int32)
        assert_array_equal(np.unique(y), [0, 1])

        model = SGDClassifier(alpha=0.1, loss="squared_hinge", n_iter=500)

        # smoke test: model is stable on scaled data
        model.fit(X_scaled, y)
        assert_true(np.isfinite(model.coef_).all())

        # model is numerically unstable on unscaled data
        msg_regxp = (
            r"Floating-point under-/overflow occurred at epoch #.*"
            " Scaling input data with StandardScaler or MinMaxScaler"
            " might help."
        )
        assert_raises_regexp(ValueError, msg_regxp, model.fit, X, y)

    def calculateLevels(self):
        """Calculate contour levels from data and settings.

        Returns levels as 1d numpy
        """

        # get dataset
        s = self.settings
        d = self.document

        minval, maxval = 0., 1.
        if s.data in d.data:
            # scan data
            data = d.data[s.data].data
            minval, maxval = N.nanmin(data), N.nanmax(data)
            if not N.isfinite(minval):
                minval = 0.
            if not N.isfinite(maxval):
                maxval = 1.

        # override if not auto
        if s.min != 'Auto':
            minval = s.min
        if s.max != 'Auto':
            maxval = s.max

        numlevels = s.numLevels
        scaling = s.scaling

        if numlevels == 1 and scaling != 'manual':
            # calculations below assume numlevels > 1
            levels = N.array([minval,])
        else:
            # trap out silly cases
            if minval == maxval:
                minval = 0.
                maxval = 1.
                
            # calculate levels for each scaling
            if scaling == 'linear':
                delta = (maxval - minval) / (numlevels-1)
                levels = minval + N.arange(numlevels)*delta
            elif scaling == 'sqrt':
                delta = N.sqrt(maxval - minval) / (numlevels-1)
                levels = minval + (N.arange(numlevels)*delta)**2
            elif scaling == 'log':
                delta = N.log(maxval/minval) / (numlevels-1)
                levels = N.exp(N.arange(numlevels)*delta)*minval
            elif scaling == 'squared':
                delta = (maxval - minval)**2 / (numlevels-1)
                levels = minval + N.sqrt(N.arange(numlevels)*delta)
            else:
                # manual
                levels = N.array(s.manualLevels)

        # for the user later
        # we do this to convert array to list of floats
        s.levelsOut = [float(i) for i in levels]

        return minval, maxval, levels

def prepare_logged(x, y):
    """
    Transform `x` and `y` to a log scale while dealing with zeros.

    This function scales `x` and `y` such that the points that are zero in one
    array are set to the min of the other array.

    When plotting expression data, frequently one sample will have reads in
    a particular feature but the other sample will not.  Expression data also
    tends to look better on a log scale, but log(0) is undefined and therefore
    cannot be shown on a plot.  This function allows these points to be shown,
    piled up along one side of the plot.

    :param x,y: NumPy arrays
    """
    xi = np.log2(x)
    yi = np.log2(y)

    xv = np.isfinite(xi)
    yv = np.isfinite(yi)

    global_min = min(xi[xv].min(), yi[yv].min())
    global_max = max(xi[xv].max(), yi[yv].max())

    xi[~xv] = global_min
    yi[~yv] = global_min

    return xi, yi

    def regularization(self, regularization):
        if regularization is None:
            self._regularization = None
            return None
        
        # Can be positive float, or positive values for all pixels.
        try:
            regularization = float(regularization)
        except (TypeError, ValueError):
            regularization = np.array(regularization).flatten()

            if regularization.size != len(self.dispersion):
                raise ValueError("regularization must be a positive value or "
                                 "an array of positive values for each pixel "
                                 "({0} != {1})".format(
                                    regularization.size,
                                    len(self.dispersion)))

            if any(0 > regularization) \
            or not np.all(np.isfinite(regularization)):
                raise ValueError("regularization terms must be "
                                 "positive and finite")
        else:
            if 0 > regularization or not np.isfinite(regularization):
                raise ValueError("regularization term must be "
                                 "positive and finite")
            regularization = np.ones_like(self.dispersion) * regularization
        self._regularization = regularization
        return None

def nextpow2(n):
    """Return the next power of 2 such as 2^p >= n.

    Notes
    -----

    Infinite and nan are left untouched, negative values are not allowed."""
    if np.any(n < 0):
        raise ValueError("n should be > 0")

    if np.isscalar(n):
        f, p = np.frexp(n)
        if f == 0.5:
            return p-1
        elif np.isfinite(f):
            return p
        else:
            return f
    else:
        f, p = np.frexp(n)
        res = f
        bet = np.isfinite(f)
        exa = (f == 0.5)
        res[bet] = p[bet]
        res[exa] = p[exa] - 1
        return res

def fft_to_hkl(h, k, l, val, coeffs, fsc_curve, resolution, full_size, flag_frac):
    '''Reformat fft record as hkl record'''
    if h or k or l:
        res = full_size / (np.linalg.norm(np.asarray([h, k, l])))
    else:
        res = 0.0

    if res < resolution or not np.isfinite(res):
        return None, None

    mag = np.abs(val)
    angle = np.angle(val, deg = True)

    if angle < 0:
        angle += 360.0

    fsc = curve_function((1. / res), coeffs, fsc_curve)
    sig = fsc_to_sigf(mag, fsc)
    fom = fsc_to_fom(fsc)
    hla, hlb = fom_to_hl(fom, np.angle(val))
    rf = bernoulli.rvs(flag_frac)
    record = np.array([h, k, l, mag, sig, angle, fom, hla, hlb, 0.0, 0.0, rf], dtype = np.float32)
    
    if not np.all(np.isfinite(record)):
        print("Skipping record %i %i %i - " %(h, k, l)),
        print(record)
        return None, None

    return record, res

def fom_to_hl(fom, phi):
    '''Convert FOMs to HLA and HLB - Kevin Cowtan - www.ysbl.york.ac.uk/~cowtan/clipper'''
    x0 = np.abs(fom)
    a0 = -7.107935 * x0
    a1 = 3.553967 - 3.524142 * x0
    a2 = 1.639294 - 2.228716 * x0
    a3 = 1.0 - x0
    w = a2 / (3.0 * a3)
    p = a1 / (3.0 * a3) - w * w
    q = -w * w * w + 0.5 * (a1 * w - a0) / a3
    d = np.sqrt(q * q + p * p * p)
    q1 = q + d
    q2 = q - d
    r1 = np.power(np.abs(q1), 1.0 / 3.0)
    r2 = np.power(np.abs(q2), 1.0 / 3.0)
    if q1 <= 0.0:
        r1 = -r1
    if q2 <= 0.0:
        r2 = -r2
    x = r1 + r2 - w
    HLA = x * np.cos(phi)
    HLB = x * np.sin(phi)
    if np.isfinite(HLA) and np.isfinite(HLB):
        return HLA, HLB
    else:
        print(" Error determining HL coefficients for FOM = "+str(fom)+' and phase = '+str(phi))
        return None, None

    def getChiImage(self, imgi=-1, img=None, srcs=None, minsb=0.):
        if img is None:
            img = self.getImage(imgi)

        # print('getChiImage:', img, ':', img.shape)
        # if srcs is None:
        #     print('Sources:')
        #     for src in self.catalog:
        #         print('  ', src)
        # else:
        #     print('Sources:', srcs)
        # print('LogPriorDerivatives:', self.getLogPriorDerivatives())
            
        mod = self.getModelImage(img, srcs=srcs, minsb=minsb)
        #print('mod:', mod.shape)
        chi = (img.getImage() - mod) * img.getInvError()
        if not np.all(np.isfinite(chi)):
            print('Chi not finite')
            print('Image finite?', np.all(np.isfinite(img.getImage())))
            print('Mod finite?', np.all(np.isfinite(mod)))
            print('InvErr finite?', np.all(np.isfinite(img.getInvError())))
            print('Current thawed parameters:')
            self.printThawedParams()
            print('Current sources:')
            for src in self.getCatalog():
                print('  ', src)
            print('Image:', img)
            print('sky:', img.getSky())
            print('psf:', img.getPsf())
        return chi

def plot(state,splits,merges):
    plt.clf()
    plt.axis([0, lifetime, -(maxl/2), maxl/2])
    for track in range(elements):
        newplot = []
        gaps=[]
        gapst=[]
        start = find_first(state[track])
        end = find_last(state[track])

        if np.isfinite(splits[track]):
            parent = splits[track]
            plt.plot([start-1,start],[state[parent][start-1][0],state[track][start][0]],'--')


        if np.isfinite(start) and np.isfinite(end):
            for t in range(start,end+1):
                newplot.append(state[track][t][0])
                if np.isfinite(state[track][t][0]): # to plot dotted lines between the gaps
                    gapst.append(t)
                    gaps.append(state[track][t][0])


            plt.plot(gapst, gaps, ':')
            plt.draw()
            plt.plot(range(start, end+1), newplot, '.-')
            plt.draw()
    plt.show()

    def __init__(self, meshsize, cellsize=None, origin=(0, 0, 0),
                 array_order=ZYX, size=None):
        if cellsize is None and size is not None:
            cellsize = np.array(
                size, dtype=float) / np.array(meshsize, dtype=float)
        self.mesh_size = np.array(meshsize, dtype=int)
        self.cell_size = np.array(cellsize, dtype=float)
        self.origin = np.array(origin, dtype=float)
        self.array_order = array_order
        self.n = np.prod(self.mesh_size)
        self.mesh_size_ao = self.mesh_size[list(array_order)]
        self.cell_size_ao = self.cell_size[list(array_order)]

        # Check validity
        assert self.mesh_size.shape == (3,)
        assert self.cell_size.shape == (3,)
        assert self.origin.shape == (3,)
        assert len(array_order) == 3

        assert self.cell_size[0] > 0 and self.cell_size[
            1] > 0 and self.cell_size[2] > 0
        assert self.mesh_size[0] > 0 and self.mesh_size[
            1] > 0 and self.mesh_size[2] > 0
        assert all(np.isfinite(self.cell_size)) and all(
            np.isfinite(self.origin))
        assert sorted(array_order) == [0, 1, 2]

 def maybe_expand_bounds(bounds):
     minval, maxval = bounds
     if not (np.isfinite(minval) and np.isfinite(maxval)):
         minval, maxval = -1.0, 1.0
     elif minval == maxval:
         minval, maxval = minval-1, minval+1
     return minval, maxval

def central_ratio(num, dnm, centerfn=np.median, finite=True):
    """Computes the central tendency (median, by default) of the ratios
    between `num` and `dnm`.  By default, this function gives the
    "Turing ratio" used in the paper by Majaj, Hong, Solomon, and DiCarlo.

    Parameters
    ----------
    num: array-like
        Numerators of ratios

    dnm: array-like, shape = `num.shape()`
        Denominators of ratios.  `num` and `dnm` must have the same shape.

    centerfn: function, optional (default=np.median)
        Function to compute the central tendency.

    finite: boolean, optional (default=True)
        If True, only finite numbers in `num` and `dnm` will be used for
        the computation of the central tendency.
    """

    num = np.array(num, dtype=DTYPE)
    dnm = np.array(dnm, dtype=DTYPE)
    assert num.shape == dnm.shape

    num = num.ravel()
    dnm = dnm.ravel()

    if finite:
        fi = np.isfinite(dnm) & np.isfinite(num)
        num = num[fi]
        dnm = dnm[fi]

    return centerfn(num / dnm)

    def start(self, f, a, b, args=()):
        r"""Prepare for the iterations."""
        self.function_calls = 0
        self.iterations = 0

        self.f = f
        self.args = args
        self.ab[:] = [a, b]
        if not np.isfinite(a) or np.imag(a) != 0:
            raise ValueError("Invalid x value: %s " % (a))
        if not np.isfinite(b) or np.imag(b) != 0:
            raise ValueError("Invalid x value: %s " % (b))

        fa = self._callf(a)
        if not np.isfinite(fa) or np.imag(fa) != 0:
            raise ValueError("Invalid function value: f(%f) -> %s " % (a, fa))
        if fa == 0:
            return _ECONVERGED, a
        fb = self._callf(b)
        if not np.isfinite(fb) or np.imag(fb) != 0:
            raise ValueError("Invalid function value: f(%f) -> %s " % (b, fb))
        if fb == 0:
            return _ECONVERGED, b

        if np.sign(fb) * np.sign(fa) > 0:
            raise ValueError("a, b must bracket a root f(%e)=%e, f(%e)=%e " %
                             (a, fa, b, fb))
        self.fab[:] = [fa, fb]

        return _EINPROGRESS, sum(self.ab) / 2.0