r"""Return the average density matrix of an inhomogeneous ensemble.
553588
def average(self, rho):
        """
        def marginal(f, rho):
            remaining = len(f.shape)
            if remaining == 0:
                return rho
            rho = sum([f[i]*rho[i] for i in range(rho.shape[0])])
            f = np.sum(f, 0)
            return marginal(f, rho)

        return marginal(self.distribution, rho)
1302506
def density(self, r, rho0, Ra, Rs):
        
        Ra, Rs = self._sort_ra_rs(Ra, Rs)
        rho = rho0 / ((1 + (r / Ra) ** 2) * (1 + (r / Rs) ** 2))
        return rho
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))
            ];
        }
    }
296299
def _sample_aAt(self,n):
        """"""
        #Sample frequency along largest eigenvalue using ARS
        dO1s=\
            bovy_ars.bovy_ars([0.,0.],[True,False],
                              [self._meandO-numpy.sqrt(self._sortedSigOEig[2]),
                               self._meandO+numpy.sqrt(self._sortedSigOEig[2])],
                              _h_ars,_hp_ars,nsamples=n,
                              hxparams=(self._meandO,self._sortedSigOEig[2]),
                              maxn=100)
        dO1s= numpy.array(dO1s)*self._sigMeanSign
        dO2s= numpy.random.normal(size=n)*numpy.sqrt(self._sortedSigOEig[1])
        dO3s= numpy.random.normal(size=n)*numpy.sqrt(self._sortedSigOEig[0])
        #Rotate into dOs in R,phi,z coordinates
        dO= numpy.vstack((dO3s,dO2s,dO1s))
        dO= numpy.dot(self._sigomatrixEig[1][:,self._sigomatrixEigsortIndx],
                      dO)
        Om= dO+numpy.tile(self._progenitor_Omega.T,(n,1)).T
        #Also generate angles
        da= numpy.random.normal(size=(3,n))*self._sigangle
        #And a random time
        dt= self.sample_t(n)
        #Integrate the orbits relative to the progenitor
        da+= dO*numpy.tile(dt,(3,1))
        angle= da+numpy.tile(self._progenitor_angle.T,(n,1)).T
        return (Om,angle,dt)
486919
def std(self, bessel_correction=True):
        """"""
        if bessel_correction:
            n = self.n
            bc = n / (n - 1)
        else:
            bc = 1
        return np.sqrt(np.average((self.bin_centers - self.mean) ** 2, weights=self.histogram)) * bc
676991
def estimate_lambda(pv):
    """"""
    LOD2 = sp.median(st.chi2.isf(pv,1))
    L = (LOD2/0.456)
    return (L)
1256683
def _num_samples(x):
    """"""
    if hasattr(x, 'fit'):
        # Don't get num_samples from an ensembles length!
        raise TypeError('Expected sequence or array-like, got '
                        'estimator %s' % x)
    if not hasattr(x, '__len__') and not hasattr(x, 'shape'):
        if hasattr(x, '__array__'):
            x = np.asarray(x)
        else:
            raise TypeError("Expected sequence or array-like, got %s" %
                            type(x))
    if hasattr(x, 'shape'):
        if len(x.shape) == 0:
            raise TypeError("Singleton array %r cannot be considered"
                            " a valid collection." % x)
        return x.shape[0]
    else:
        return len(x)
381957
def mag(z):
    """"""
    if isinstance(z[0], np.ndarray):
        return np.array(list(map(np.linalg.norm, z)))
    else:
        return np.linalg.norm(z)
1188507
def pvalues(self):
        """"""
        self.compute_statistics()

        lml_alts = self.alt_lmls()
        lml_null = self.null_lml()

        lrs = -2 * lml_null + 2 * asarray(lml_alts)

        from scipy.stats import chi2
        chi2 = chi2(df=1)

        return chi2.sf(lrs)
797196
def _density_par(self,dangle,tdisrupt=None):
        """"""
        if tdisrupt is None: tdisrupt= self._tdisrupt
        dOmin= dangle/tdisrupt
        # Normalize to 1 close to progenitor
        return 0.5\
            *(1.+special.erf((self._meandO-dOmin)\
                                 /numpy.sqrt(2.*self._sortedSigOEig[2])))
194492
def _pdf(self, x, dist, cache):
        """"""
        return evaluation.evaluate_density(
            dist, numpy.arcsinh(x), cache=cache)/numpy.sqrt(1+x*x)
213397
def _get_kriging_matrix(self, n, n_withdrifts):
        """"""

        xyz = np.concatenate((self.X_ADJUSTED[:, np.newaxis], self.Y_ADJUSTED[:, np.newaxis],
                              self.Z_ADJUSTED[:, np.newaxis]), axis=1)
        d = cdist(xyz, xyz, 'euclidean')
        if self.UNBIAS:
            a = np.zeros((n_withdrifts+1, n_withdrifts+1))
        else:
            a = np.zeros((n_withdrifts, n_withdrifts))
        a[:n, :n] = - self.variogram_function(self.variogram_model_parameters, d)
        np.fill_diagonal(a, 0.)

        i = n
        if self.regional_linear_drift:
            a[:n, i] = self.X_ADJUSTED
            a[i, :n] = self.X_ADJUSTED
            i += 1
            a[:n, i] = self.Y_ADJUSTED
            a[i, :n] = self.Y_ADJUSTED
            i += 1
            a[:n, i] = self.Z_ADJUSTED
            a[i, :n] = self.Z_ADJUSTED
            i += 1
        if self.specified_drift:
            for arr in self.specified_drift_data_arrays:
                a[:n, i] = arr
                a[i, :n] = arr
                i += 1
        if self.functional_drift:
            for func in self.functional_drift_terms:
                a[:n, i] = func(self.X_ADJUSTED, self.Y_ADJUSTED, self.Z_ADJUSTED)
                a[i, :n] = func(self.X_ADJUSTED, self.Y_ADJUSTED, self.Z_ADJUSTED)
                i += 1
        if i != n_withdrifts:
            warnings.warn("Error in creating kriging matrix. Kriging may fail.", RuntimeWarning)
        if self.UNBIAS:
            a[n_withdrifts, :n] = 1.0
            a[:n, n_withdrifts] = 1.0
            a[n:n_withdrifts + 1, n:n_withdrifts + 1] = 0.0

        return a
394987
def _M(self, k, t, tips=None, gaps=None):
        """"""
        if self._distributionmodel:
            return self.model.M(k, t, tips, gaps)
        else:
            return self.model.M(t, tips, gaps)
623725
def matrix(self, _):
        """"""
        x, y, z = self.coord
        return np.array([
            [1., 0., 0., x],
            [0., 1., 0., y],
            [0., 0., 1., z],
            [0., 0., 0., 1.]
        ])
723993
def density_2d(self, x, y, rho0, gamma, center_x=0, center_y=0):
        
        x_ = x - center_x
        y_ = y - center_y
        r = np.sqrt(x_**2 + y_**2)
        sigma = np.sqrt(np.pi) * special.gamma(1./2*(-1+gamma))/special.gamma(gamma/2.) * r**(1-gamma) * rho0
        return sigma
731013
def getCovariance(self,normalize=True,i0=None,i1=None,pos0=None,pos1=None,chrom=None,center=True,unit=True,pos_cum0=None,pos_cum1=None,blocksize=None,X=None,**kw_args):
        """"""
        if X is not None:
            K=X.dot(X.T)
            Nsnp=X.shape[1]
        else:
            if (i0 is None) and (i1 is None) and ((pos0 is not None) & (pos1 is not None) & (chrom is not None)) or ((pos_cum0 is not None) & (pos_cum1 is not None)):
                i0,i1=self.getGenoIndex(pos0=pos0,pos1=pos1,chrom=chrom,pos_cum0=pos_cum0,pos_cum1=pose_cum1)

            [N,M]=self.genoM.shape
            if blocksize is None:
                blocksize=M
            if i0 is None:
                i0=0
            if i1 is None:
                i1=M
            nread = i0
            K=None
            Nsnp=i1-i0
            while nread<i1:
                thisblock=min(blocksize,i1-nread)
                X=self.getGenotypes(i0=nread,i1=(nread+thisblock),center=center,unit=unit,**kw_args)    
                if K is None:
                    K=X.dot(X.T)
                else:
                    K+=X.dot(X.T)
                nread+=thisblock
        if normalize:
            K/=(K.diagonal().mean())
        else:#divide by number of SNPs in K
            K/=Nsnp
        return K
71573
def _density_par_approx_higherorder(self,Oparb,lowbindx,
                                        _return_array=False,
                                        gaussxpolyInt=None):
        """"""
        spline_order= self._kick_interpdOpar_raw._eval_args[2]
        if spline_order == 1: return 0.
        # Form all Gaussian-like integrals necessary
        ll= (numpy.roll(Oparb,-1)[:-1]-self._kick_interpdOpar_poly.c[-1]\
            -self._meandO\
            -self._kick_interpdOpar_poly.c[-2]*self._timpact\
            *(Oparb-numpy.roll(Oparb,-1))[:-1])\
            /numpy.sqrt(2.*self._sortedSigOEig[2])
        ul= (Oparb[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO)\
            /numpy.sqrt(2.*self._sortedSigOEig[2])
        if gaussxpolyInt is None:
            gaussxpolyInt=\
                self._densMoments_approx_higherorder_gaussxpolyInts(\
                ll,ul,spline_order+1)
        # Now multiply in the coefficients for each order
        powers= numpy.tile(numpy.arange(spline_order+1)[::-1],
                           (len(ul),1)).T
        gaussxpolyInt*= -0.5*(-numpy.sqrt(2.))**(powers+1)\
            *self._sortedSigOEig[2]**(0.5*(powers-1))
        powers= numpy.tile(numpy.arange(spline_order+1)[::-1][:-2],
                           (len(ul),1)).T
        for jj in range(spline_order+1):
            gaussxpolyInt[-jj-1]*= numpy.sum(\
                self._kick_interpdOpar_poly.c[:-2]
                *self._timpact**powers
                /(1.+self._kick_interpdOpar_poly.c[-2]*self._timpact)
                **(powers+1)
                *special.binom(powers,jj)
                *(Oparb[:-1]-self._kick_interpdOpar_poly.c[-1]-self._meandO)
                **(powers-jj),axis=0)
        if _return_array:
            return numpy.sum(gaussxpolyInt,axis=0)
        else:
            return numpy.sum(gaussxpolyInt[:,:lowbindx+1])