Kolmorgorov-Smirnov statistic. Finds the
    probability that the data are distributed
    as func - used method of Numerical Recipes (Press et al., 1986)
778417
def k_s(X):
    
    xbar, sigma = pmag.gausspars(X)
    d, f = 0, 0.
    for i in range(1, len(X) + 1):
        b = old_div(float(i), float(len(X)))
        a = gaussfunc(X[i - 1], xbar, sigma)
        if abs(f - a) > abs(b - a):
            delta = abs(f - a)
        else:
            delta = abs(b - a)
        if delta > d:
            d = delta
        f = b
    return d, xbar, sigma
327105
def _statsmodels_bivariate_kde(self, x, y, bws, xsize, ysize, xyranges):
        """"""
        import statsmodels.nonparametric.api as smnp
        for i, (coord, bw) in enumerate(zip([x, y], bws)):
            if isinstance(bw, six.string_types):
                bw_func = getattr(smnp.bandwidths, "bw_" + bw)
                bws[i] = bw_func(coord)
        kde = smnp.KDEMultivariate([x, y], "cc", bws)
        x_support = np.linspace(xyranges[0][0], xyranges[0][1], xsize)
        y_support = np.linspace(xyranges[1][0], xyranges[1][1], ysize)
        xx, yy = np.meshgrid(x_support, y_support)
        z = kde.pdf([xx.ravel(), yy.ravel()]).reshape(xx.shape)
        return x_support, y_support, z
455576
function gaussianKernel1D(s, n) {
    if (typeof n !== 'number') n = 7;
    var kernel = new Float64Array(n),
        halfN = Math.floor(n * 0.5),
        odd = n % 2,
        i;
    if (!s || !n) return kernel;
    for (i = 0; i <= halfN; i++) {
        kernel[i] = gauss(s * (i - halfN - odd * 0.5), s);
    }
    for (; i < n; i++) {
        kernel[i] = kernel[n - 1 - i];
    }
    return kernel;
}
883182
def sample_vMF(mu, kappa, num_samples):
    
    dim = len(mu)
    result = np.zeros((num_samples, dim))
    for nn in range(num_samples):
        # sample offset from center (on sphere) with spread kappa
        w = _sample_weight(kappa, dim)

        # sample a point v on the unit sphere that's orthogonal to mu
        v = _sample_orthonormal_to(mu)

        # compute new point
        result[nn, :] = v * np.sqrt(1. - w ** 2) + w * mu

    return result
156222
def bootstrap_noise(data, func, n=10000, std=1, symmetric=True):
    
    boot_dist = []
    arr = N.zeros(data.shape)
    for i in range(n):
        if symmetric:
            # Noise on all three axes
            arr = N.random.randn(*data.shape)*std
        else:
            # Only z-coordinate noise
            arr[:,-1] = N.random.randn(data.shape[0])*std
        boot_dist.append(func(data+arr))
    return N.array(boot_dist)
867241
def _dEndfr(self):
        
        Mc = self._chirp_mass()
        return (np.pi**(2./3.)*Mc**(5./3.)/(3.*(1.+self.z)**(1./3.)
                * (self.freqs_orb/(1.+self.z))**(1./3.))*(2./self.n)**(2./3.)
                * self._g_func()/self._f_func())
772271
def cross_correlate(observed, template_dispersion, template_fluxes,
    rebin="template", wavelength_range=None, continuum_degree=-1,
    apodize=0.10, rescale=False, z_limits=None, full_output=False):
    

    template_fluxes = np.atleast_2d(template_fluxes)
    template_dispersion = np.array(template_dispersion)
    if template_dispersion.shape[0] != template_fluxes.shape[1]:
        raise ValueError("template dispersion must have size (N_pixels,) and "\
            "template fluxes must have size (N_models, N_pixels)")

    try:
        continuum_degree = int(continuum_degree)
    except (TypeError, ValueError):
        raise TypeError("continuum order must be an integer-like object")

    if not (1 > apodize >= 0):
        raise ValueError("apodize fraction must be between 0 and 1")

    if rebin.lower() == "template":
        # Put the template fluxes onto the observed dispersion map.
        dispersion = observed.disp
        template_flux = template_fluxes \
            * resample(template_dispersion, observed.disp)
        observed_flux = observed.flux.copy()
        observed_ivar = observed.ivariance.copy()

    elif rebin.lower() == "observed":
        # Put the observed fluxes onto the template dispersion map.
        dispersion = template_dispersion
        template_flux = template_fluxes
        mat = resample(observed.disp, template_dispersion)
        observed_flux = observed.flux.copy() * mat
        observed_ivar = observed.ivariance.copy() * mat
        
    else:
        raise ValueError("rebin must be either `template` or `observed`")

    if wavelength_range is not None:
        if not isinstance(wavelength_range, (tuple, list, np.ndarray)) \
        or len(wavelength_range) != 2:
            raise TypeError("wavelength range must be a two length tuple")

        indices = np.clip(
            template_dispersion.searchsorted(wavelength_range) + [0, 1],
            0, template_dispersion.size)
        template_dispersion = template_dispersion[indices[0]:indices[1]]
        template_flux = template_flux[indices[0]:indices[1]]
        observed_flux = observed_flux[indices[0]:indices[1]]

    N_templates = template_fluxes.shape[0]

    # Clip out the non-finite edges, if they exist.
    _ = np.where(np.isfinite(template_flux[0] * observed_flux))[0]
    l_idx, u_idx = _[0], _[-1]

    dispersion = dispersion[l_idx:u_idx]
    observed_flux = observed_flux[l_idx:u_idx]
    observed_ivar = observed_ivar[l_idx:u_idx]
    template_flux = template_flux[:, l_idx:u_idx]

    # Ensure an even number of points.
    N = u_idx - l_idx
    N = N - 1 if N % 2 > 0 else N
    dispersion = dispersion[:N]
    observed_flux = observed_flux[:N]
    observed_ivar = observed_ivar[:N]
    template_flux = template_flux[:, :N]

    finite = np.isfinite(observed_flux * observed_ivar) * (observed_flux > 1e-3)
    
    # Continuum.
    if continuum_degree >= 0:
        coefficients = np.polyfit(dispersion[finite], observed_flux[finite],
            continuum_degree)
        observed_flux /= np.polyval(coefficients, dispersion)

    # Interpolate over non-finite pixels.
    observed_flux[~finite] = np.interp(dispersion[~finite], dispersion[finite],
        observed_flux[finite], left=1, right=1)

    # Scale the flux level to that the template intensities
    if rescale:
        observed_flux = (observed_flux * template_flux.ptp()) \
            + template_flux.min()

    # Apodize edges.
    edge_buffer = apodize * (dispersion[-1] - dispersion[0])
    low_w_indices = np.nonzero(dispersion < dispersion[0] + edge_buffer)[0]
    high_w_indices = np.nonzero(dispersion > dispersion[-1] - edge_buffer)[0]

    apod_curve = np.ones(N, dtype='d')
    apod_curve[low_w_indices] = (1.0 + np.cos(np.pi*(
        1.0 - (dispersion[low_w_indices] - dispersion[0])/edge_buffer)))/2.
    apod_curve[high_w_indices] = (1.0 + np.cos(np.pi*(
        1.0 - (dispersion[-1] - dispersion[high_w_indices])/edge_buffer)))/2.

    apod_observed_flux = observed_flux * apod_curve
    apod_template_flux = template_flux * apod_curve

    fft_observed_flux = np.fft.fft(apod_observed_flux)
    fft_template_flux = np.fft.fft(apod_template_flux)
    template_flux_corr = (fft_observed_flux * fft_template_flux.conjugate())
    template_flux_corr /= np.sqrt(
        np.inner(apod_observed_flux, apod_observed_flux))

    z_array = np.array(dispersion.copy())/dispersion[N/2] - 1.0

    z = np.ones(N_templates) * np.nan
    z_err = np.ones(N_templates) * np.nan
    R = np.ones(N_templates) * np.nan
    for i in range(N_templates):

        denominator = np.sqrt(
            np.inner(apod_template_flux[i, :], apod_template_flux[i, :]))
        flux_correlation = template_flux_corr[i, :] / denominator 
        correlation = np.fft.ifft(flux_correlation).real

        # Reflect about zero
        ccf = np.zeros(N)
        ccf[:N/2] = correlation[N/2:]
        ccf[N/2:] = correlation[:N/2]

        # Get height and redshift of best peak
        h = ccf.max()

        # Scale the CCF
        #ccf -= ccf.min()
        #ccf *= (h/ccf.max())
        
        R[i] = h
        z[i] = z_array[ccf.argmax()]
        try:
            z_err[i] = (np.ptp(z_array[np.where(ccf >= 0.5*h)])/2.35482)**2
        except ValueError:
            continue

    # Limits:
    if z_limits is not None:
        R[~((z_limits[1] >= z) * (z >= z_limits[0]))] = np.nan

    # Re-measure the velocity at the best peak.
    index = np.nanargmax(R)

    denominator = np.sqrt(
        np.inner(apod_template_flux[index, :], apod_template_flux[index, :]))

    flux_correlation = template_flux_corr[index, :] / denominator
    correlation = np.fft.ifft(flux_correlation).real

    # FReflect about zero.
    ccf = np.zeros(N)
    ccf[:N/2] = correlation[N/2:]
    ccf[N/2:] = correlation[:N/2]

    h = ccf.max()
    ccf -= ccf.min()
    ccf *= h/ccf.max()

    # Fit +/- 5 pixels
    idx = np.argmax(ccf) - 3, np.argmax(ccf) + 3

    coeffs = np.polyfit(z_array[idx[0]:idx[1]], ccf[idx[0]:idx[1]], 2)
    x_i = np.linspace(z_array[idx[0]], z_array[idx[1]], 1000)
    y_i = np.polyval(coeffs, x_i)

    # Update the value
    z[index] = x_i[y_i.argmax()]

    """
    # Fit the profile peak.
    fwhm = np.ptp(z_array[np.where(ccf >= 0.5)])
    p0 = np.array([z[index], fwhm/2.355, 1.0]) # mu, sigma, peak.
    use = (z_array > (p0[0] - 3 * p0[1])) * ((p0[0] + 3*p0[1]) > z_array)
    x, y = z_array[use], ccf[use]

    # Fit profile:
    f = lambda x, mu, sigma, peak: peak * np.exp(-(mu - x)**2/(2*sigma**2)) 
    
    import matplotlib.pyplot as plt
    plt.plot(z_array, ccf, c='k')
    plt.plot(x, y, c='b')
    plt.plot(x, f(x, *p0), c='r')

    p1 = curve_fit(f, x, y, p0=p0)
    plt.plot(x, f(x, *p1[0]), c='m')

    plt.gca().set_ylim(0,1)
    raise a

    p1 = curve_fit(f, x, y, p0=p0)
    """

    return (z * c, z_err * c, R)
103843
function (data) {
    var histogram = data || this.get('data');
    var freqAccessor = function (a) { return a.freq; };
    var osc = d3.max(histogram, freqAccessor) - d3.min(histogram, freqAccessor);
    var mean = d3.mean(histogram, freqAccessor);
    // When the difference between the max and the min values is less than
    // 10 percent of the mean, it's a flat histogram (F)
    if (osc < mean * 0.1) return 'F';
    var sumFreqs = d3.sum(histogram, freqAccessor);
    var freqs = histogram.map(function (bin) {
      return 100 * bin.freq / sumFreqs;
    });

    // The ajus array represents relative growths
    var ajus = freqs.map(function (freq, index) {
      var next = freqs[index + 1];
      if (freq > next) return -1;
      if (Math.abs(freq - next) <= 0.05) return 0;
      return 1;
    });
    ajus.pop();
    var maxAjus = d3.max(ajus);
    var minAjus = d3.min(ajus);
    // If it never grows or shrinks, it returns flat
    if (minAjus === 0 && maxAjus === 0) return 'F';
    else if (maxAjus < 1) return 'L';
    else if (minAjus > -1) return 'J';
    else {
      var uniques = _.uniq(ajus);
      var A_TYPES = [[1, -1], [1, 0, -1], [1, -1, 0], [0, 1, -1]];
      var U_TYPES = [[-1, 1], [-1, 0, 1], [-1, 1, 0], [0, -1, 1]];
      if (A_TYPES.some(function (e) {
        return _.isEqual(e, uniques);
      })) return 'A';
      else if (U_TYPES.some(function (e) {
        return _.isEqual(e, uniques);
      })) return 'U';
      else return 'S';
    }
  }
1064109
def show_G_distribution(data):
    ''''''
    Xs, t = fitting.preprocess_data(data)  

    Theta, Phi = np.meshgrid(np.linspace(0, np.pi, 50), np.linspace(0, 2 * np.pi, 50))
    G = []

    for i in range(len(Theta)):
        G.append([])
        for j in range(len(Theta[i])):
            w = fitting.direction(Theta[i][j], Phi[i][j])
            G[-1].append(fitting.G(w, Xs))

    plt.imshow(G, extent=[0, np.pi, 0, 2 * np.pi], origin='lower')
    plt.show()
880466
def poisson_consensus_se(data, k, n_runs=10, **se_params):
    
    clusters = []
    for i in range(n_runs):
        assignments, means  = poisson_cluster(data, k)
        clusters.append(assignments)
    clusterings = np.vstack(clusters)
    consensus = CE.cluster_ensembles(clusterings, verbose=False, N_clusters_max=k)
    init_m, init_w = nmf_init(data, consensus, k, 'basic')
    M, W, ll = poisson_estimate_state(data, k, init_means=init_m, init_weights=init_w, **se_params)
    return M, W, ll
956899
function(distribution, coverageAlpha) {
        if (distribution.numTrials < 1000) {
            // don't even bother trying to optimize for small-ish sample sizes
            return [0, distribution.numTrials];
        } else {
            return [
                Math.floor(distribution.inverseCdf(coverageAlpha / 2)),
                Math.ceil(distribution.inverseSurvival(coverageAlpha / 2))
            ];
        }
    }
93969
def bingham_pdf(fit):
    
    # Uses eigenvectors of the covariance matrix
    e = fit.hyperbolic_axes #singular_values
    #e = sampling_covariance(fit) # not sure
    e = e[2]**2/e

    kappa = (e-e[2])[:-1]
    kappa /= kappa[-1]
    F = N.sqrt(N.pi)*confluent_hypergeometric_function(*kappa)


    ax = fit.axes

    Z = 1/e
    M = ax
    F = 1/hyp1f1(*1/Z)

    def pdf(coords):
        lon,lat = coords

        I = lat
        D = lon# + N.pi/2
        #D,I = _rotate(N.degrees(D),N.degrees(I),90)

        # Bingham is given in spherical coordinates of inclination
        # and declination in radians

        # From USGS bingham statistics reference

        xhat = N.array(sph2cart(lon,lat)).T

        #return F*expm(dot(xhat.T, M, N.diag(Z), M.T, xhat))

        return 1/(F*N.exp(dot(xhat.T, M, N.diag(Z), M.T, xhat)))

    return pdf
1028992
def numbaGaussian2d(psf, sy, sx):
    '''
    
    '''
    ps0, ps1 = psf.shape
    c0,c1 = ps0//2, ps1//2
    ssx = 2*sx**2
    ssy = 2*sy**2
    for i in range(ps0):
        for j in range(ps1):
            psf[i,j]=exp( -( (i-c0)**2/ssy
                            +(j-c1)**2/ssx) )
    psf/=psf.sum()
669550
def _get_stddevs(self, C, stddev_types, num_sites):
        
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        if self.DEFINED_FOR_INTENSITY_MEASURE_COMPONENT == 'Random horizontal':
            # Using equation 8 of the paper,
            # corrected as indicated in the erratum
            Sr = np.sqrt(C['SlZ']**2 + (C['S3']/np.sqrt(2))**2)
        else:
            Sr = C['SlZ']

        stddevs = [np.zeros(num_sites) + Sr for _ in stddev_types]
        return stddevs
108202
def obfn_dfd(self):
        
        """

        Ef = sl.inner(self.Zf, self.obfn_fvarf(), axis=self.cri.axisM) \
          - self.Sf
        return (np.linalg.norm(self.W * sl.irfftn(Ef, self.cri.Nv,
                                                  self.cri.axisN))**2) / 2.0
406880
def calcSMAfromT(self, epsilon=0.7):
        

        return eq.MeanPlanetTemp(self.albedo(), self.star.T, self.star.R, epsilon, self.T).a
934869
def hessian(self, x, y, Rs, theta_Rs, r_core, center_x=0, center_y=0):

        

        if Rs < 0.0001:
            Rs = 0.0001
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_ ** 2 + y_ ** 2)

        rho0 = self._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs, r_core=r_core)

        kappa = self.density_2d(x_, y_, Rs, rho0, r_core)

        gamma1, gamma2 = self.cBurkGamma(R, Rs, rho0, r_core, x_, y_)
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy