Python scipy 模块，pi() 实例源码

def rvs( self ):
        if self.ss == 0:
            x = random.random() * self.w
            y = random.random() * self.h
            self.set_point( x, y )
        while self.qs:
            x_idx = int( random.random() * self.qs )
            s = self.queue[ x_idx ]
            for y_idx in range( self.n ):
                a = 2 * scipy.pi * random.random()
                b = scipy.sqrt( self.A * random.random() + self.r2 )
                x = s[0] + b*scipy.cos( a )
                y = s[1] + b*scipy.sin( a )
                if( x >= 0 )and( x < self.w ):
                    if( y >= 0 )and( y < self.h ):
                        if( self.distance( x, y ) ):
                            self.set_point( x, y )
            del self.queue[x_idx]
            self.qs -= 1
        sample = list( filter( None, self.grid ) )
        sample = scipy.asfarray( sample )
        return sample

def fit(self, X, w):
        if len(X) == 0:
            raise NotEnoughParticles("Fitting not possible.")
        self.X_arr = X.as_matrix()

        ctree = cKDTree(X)
        _, indices = ctree.query(X, k=min(self.k + 1, X.shape[0]))

        covs, inv_covs, dets = list(zip(*[self._cov_and_inv(n, indices)
                                    for n in range(X.shape[0])]))
        self.covs = sp.array(covs)
        self.inv_covs = sp.array(inv_covs)
        self.determinants = sp.array(dets)

        self.normalization = sp.sqrt(
            (2 * sp.pi) ** self.X_arr.shape[1] * self.determinants)

        if not sp.isreal(self.normalization).all():
            raise Exception("Normalization not real")
        self.normalization = sp.real(self.normalization)

def model(THETA, time, kplanets):
    modelo = 0.0
    if kplanets == 0:
        return 0.0
    for i in range(kplanets):
        As, P, Ac, S, C = THETA[5*i:5*(i+1)]
        A = As ** 2 + Ac ** 2
        ecc = S ** 2 + C ** 2
        w = sp.arccos(C / (ecc ** 0.5))  # longitude of periastron
        phase = sp.arccos(Ac / (A ** 0.5))
        ### test
        if S < 0:
            w = 2 * sp.pi - sp.arccos(C / (ecc ** 0.5))
        if As < 0:
            phase = 2 * sp.pi - sp.arccos(Ac / (A ** 0.5))
        ###
        per = sp.exp(P)
        freq = 2. * sp.pi / per
        M = freq * time + phase
        E = sp.array([MarkleyKESolver().getE(m, ecc) for m in M])
        f = (sp.arctan(((1. + ecc) ** 0.5 / (1. - ecc) ** 0.5) * sp.tan(E / 2.)) * 2.)
        modelo += A * (sp.cos(f + w) + ecc * sp.cos(w))
    return  modelo

def kde(data, N=None, MIN=None, MAX=None):

    # Parameters to set up the mesh on which to calculate
    N = 2**14 if N is None else int(2**sci.ceil(sci.log2(N)))
    if MIN is None or MAX is None:
        minimum = min(data)
        maximum = max(data)
        Range = maximum - minimum
        MIN = minimum - Range/10 if MIN is None else MIN
        MAX = maximum + Range/10 if MAX is None else MAX

    # Range of the data
    R = MAX-MIN

    # Histogram the data to get a crude first approximation of the density
    M = len(data)
    DataHist, bins = sci.histogram(data, bins=N, range=(MIN,MAX))
    DataHist = DataHist/M
    DCTData = scipy.fftpack.dct(DataHist, norm=None)

    I = [iN*iN for iN in xrange(1, N)]
    SqDCTData = (DCTData[1:]/2)**2

    # The fixed point calculation finds the bandwidth = t_star
    guess = 0.1
    try:
        t_star = scipy.optimize.brentq(fixed_point, 0, guess, 
                                       args=(M, I, SqDCTData))
    except ValueError:
        print 'Oops!'
        return None

    # Smooth the DCTransformed data using t_star
    SmDCTData = DCTData*sci.exp(-sci.arange(N)**2*sci.pi**2*t_star/2)
    # Inverse DCT to get density
    density = scipy.fftpack.idct(SmDCTData, norm=None)*N/R
    mesh = [(bins[i]+bins[i+1])/2 for i in xrange(N)]
    bandwidth = sci.sqrt(t_star)*R

    density = density/sci.trapz(density, mesh)
    return bandwidth, mesh, density

def fixed_point(t, M, I, a2):
    l=7
    I = sci.float128(I)
    M = sci.float128(M)
    a2 = sci.float128(a2)
    f = 2*sci.pi**(2*l)*sci.sum(I**l*a2*sci.exp(-I*sci.pi**2*t))
    for s in range(l, 1, -1):
        K0 = sci.prod(xrange(1, 2*s, 2))/sci.sqrt(2*sci.pi)
        const = (1 + (1/2)**(s + 1/2))/3
        time=(2*const*K0/M/f)**(2/(3+2*s))
        f=2*sci.pi**(2*s)*sci.sum(I**s*a2*sci.exp(-I*sci.pi**2*time))
    return t-(2*M*sci.sqrt(sci.pi)*f)**(-2/5)

def nLLeval(ldelta, Uy, S, REML=True):
    """
    evaluate the negative log likelihood of a random effects model:
    nLL = 1/2(n_s*log(2pi) + logdet(K) + 1/ss * y^T(K + deltaI)^{-1}y,
    where K = USU^T.

    Uy: transformed outcome: n_s x 1
    S:  eigenvectors of K: n_s
    ldelta: log-transformed ratio sigma_gg/sigma_ee
    """
    n_s = Uy.shape[0]
    delta = scipy.exp(ldelta)

    # evaluate log determinant
    Sd = S + delta
    ldet = scipy.sum(scipy.log(Sd))

    # evaluate the variance
    Sdi = 1.0 / Sd
    Sdi=Sdi.reshape((Sdi.shape[0],1))
    ss = 1. / n_s * (Uy*Uy*Sdi).sum()
    # evalue the negative log likelihood
    nLL = 0.5 * (n_s * scipy.log(2.0 * scipy.pi) + ldet + n_s + n_s * scipy.log(ss))

    if REML:
        pass

    return nLL

def ellipse2bbox(a, b, angle, cx, cy):
        a, b = max(a, b), min(a, b)
        ca = sp.cos(angle)
        sa = sp.sin(angle)
        if sa == 0.0:
            cta = 2.0 / sp.pi
        else:
            cta = ca / sa

        if ca == 0.0:
            ta = sp.pi / 2.0
        else:
            ta = sa / ca

        x = lambda t: cx + a * sp.cos(t) * ca - b * sp.sin(t) * sa


        y = lambda t: cy + b * sp.sin(t) * ca + a * sp.cos(t) * sa

        # x = cx + a * cos(t) * cos(angle) - b * sin(t) * sin(angle)
        # tan(t) = -b * tan(angle) / a
        tx1 = sp.arctan(-b * ta / a)
        tx2 = tx1 - sp.pi
        x1, y1 = x(tx1), y(tx1)
        x2, y2 = x(tx2), y(tx2)

        # y = cy + b * sin(t) * cos(angle) + a * cos(t) * sin(angle)
        # tan(t) = b * cot(angle) / a
        ty1 = sp.arctan(b * cta / a)
        ty2 = ty1 - sp.pi
        x3, y3 = x(ty1), y(ty1)
        x4, y4 = x(ty2), y(ty2)

        minx, maxx = Util.minmax([x1, x2, x3, x4])
        miny, maxy = Util.minmax([y1, y2, y3, y4])
        return sp.floor(minx), sp.floor(miny), sp.ceil(maxx), sp.ceil(maxy)

def _pipePositions(Ds, nPipes):
    """ Positions pipes in an axisymetric configuration.
    """
    dt = pi / float(nPipes)
    pos = [(0., 0.) for i in range(2*nPipes)]
    for i in range(nPipes):
        pos[i] = (Ds*np.cos(2.0*i*dt+pi), Ds*np.sin(2.0*i*dt+pi))
        pos[i+nPipes] = (Ds*np.cos(2.0*i*dt+pi+dt), Ds*np.sin(2.0*i*dt+pi+dt))
    return pos


# Main function

def semimodel(self, params, time):
        A, P, phase, w, ecc = params
        freq = 2. * sp.pi / P
        M = freq * time + phase
        E = sp.array([MarkleyKESolver().getE(m, ecc) for m in M])
        f = (sp.arctan(((1. + ecc) ** 0.5 / (1. - ecc) ** 0.5) * sp.tan(E / 2.)) * 2.)
        modelo = A * (sp.cos(f + w) + ecc * sp.cos(w))
        return modelo

def henshin(self, thetas, kplanets):
        for i in range(kplanets):
            Ask = thetas[:, i*5]
            Pk = thetas[:, i*5 + 1]
            Ack = thetas[:, i*5 + 2]
            Sk = thetas[:, i*5 + 3]
            Ck = thetas[:, i*5 + 4]

            ecck  = Sk ** 2 + Ck ** 2
            Ak = Ask ** 2 + Ack ** 2
            wk = sp.arccos(Ck / (ecck ** 0.5))
            Phasek = sp.arccos(Ack / (Ak ** 0.5))
            for j in range(len(Sk)):
                if Sk[j] < 0:
                    wk[j] = 2 * sp.pi - sp.arccos(Ck[j] / (ecck[j] ** 0.5))
                if Ask[j] < 0:
                    Phasek[j] = 2 * sp.pi - sp.arccos(Ack[j] / (Ak[j] ** 0.5))

            thetas[:, i*5] = Ak
            thetas[:, i*5 + 1] = sp.exp(Pk)
            thetas[:, i*5 + 2] = Phasek
            thetas[:, i*5 + 3] = wk
            thetas[:, i*5 + 4] = ecck
        return thetas

    ########################################
    # los IC sirven para los mod_lims? NO!
    # mod_lims sólo acotan el espacio donde buscar, excepto para periodo

def nLLeval(ldelta, Uy, S, REML=True):
    """
    evaluate the negative log likelihood of a random effects model:
    nLL = 1/2(n_s*log(2pi) + logdet(K) + 1/ss * y^T(K + deltaI)^{-1}y,
    where K = USU^T.
    Uy: transformed outcome: n_s x 1
    S:  eigenvectors of K: n_s
    ldelta: log-transformed ratio sigma_gg/sigma_ee
    """
    n_s = Uy.shape[0]
    delta = scipy.exp(ldelta)

    # evaluate log determinant
    Sd = S + delta
    ldet = scipy.sum(scipy.log(Sd))

    # evaluate the variance
    Sdi = 1.0 / Sd
    Uy = Uy.flatten()
    ss = 1. / n_s * (Uy * Uy * Sdi).sum()

    # evalue the negative log likelihood
    nLL = 0.5 * (n_s * scipy.log(2.0 * scipy.pi) + ldet + n_s + n_s * scipy.log(ss))

    if REML:
        pass

    return nLL

def func_2vN(Ek, Ek_grid, l, eta, hfk):
    """
    Linearly interpolate the value of hfk on Ek_grid at point Ek.

    Parameters
    ----------
    Ek : float
        Energy value (not necessarily a grid point).
    Ek_grid : array
        Energy grid.
    l : int
        Lead label.
    eta : int
        Integer describing (+/-1) if infinitesimal eta is positive or negative.
    hfk : array
        Array containing Hilbert transform of Fermi function (or 1-Fermi).

    Returns
    -------
    float
        Interpolated value of hfk at Ek.
    """
    if Ek<Ek_grid[0] or Ek>Ek_grid[-1]:
        return 0
    #
    b_idx = int((Ek-Ek_grid[0])/(Ek_grid[1]-Ek_grid[0]))+1
    if b_idx == len(Ek_grid): b_idx -= 1
    a_idx = b_idx - 1
    b, a = Ek_grid[b_idx], Ek_grid[a_idx]
    #
    fb = hfk[l, b_idx]
    fa = hfk[l, a_idx]
    rez = (fb-fa)/(b-a)*Ek + (b*fa-a*fb)/(b-a)
    return pi*rez if eta+1 else pi*rez.conjugate()

def func_pauli(Ecb, mu, T, Dm, Dp, itype):
    """
    Function used when generating Pauli master equation kernel.

    Parameters
    ----------
    Ecb : float
        Energy.
    mu : float
        Chemical potential.
    T : float
        Temperature.
    Dm,Dp : float
        Bandwidth.

    Returns
    -------
    array
        | Array of two float numbers [cur0, cur1] containing
          momentum-integrated current amplitudes.
        | cur0 - particle current amplitude.
        | cur1 - hole current amplitude.
    """
    alpha = (Ecb-mu)/T
    Rm, Rp = (Dm-mu)/T, (Dp-mu)/T
    if itype == 1 or itype == 3 or (alpha < Rp and alpha > Rm):
        cur0 = fermi_func(alpha)
        cur1 = 1-cur0
        rez = 2*pi*np.array([cur0, cur1])
    else:
        rez = np.zeros(2)
    return rez

def logl(theta, time, rv, err, ins, staract, starflag, kplanets, nins, MOAV, totcornum):
    i, lnl = 0, 0
    ndat = len(time)
    model_params = kplanets * 5
    ins_params = nins * 2 * (MOAV + 1)
    jitter, offset, macoef, timescale = sp.zeros(ndat), sp.zeros(ndat), sp.array([sp.zeros(ndat) for i in range(MOAV)]), sp.array([sp.zeros(ndat) for i in range(MOAV)])
    ACC = theta[model_params] * (time - time[0])

    residuals = sp.zeros(ndat)
    for i in range(ndat):
        jitpos = int(model_params + ins[i] * 2 * (MOAV+1) + 1)
        jitter[i], offset[i] = theta[jitpos], theta[jitpos + 1]  # jitt
        for j in range(MOAV):
            macoef[j][i], timescale[j][i] = theta[jitpos + 2*(j+1)], theta[jitpos + 2*(j+1) + 1]
    a1 = (theta[:model_params])

    if totcornum:
        COR = sp.array([sp.array([sp.zeros(ndat) for k in range(len(starflag[i]))]) for i in range(len(starflag))])
        SA = theta[model_params+ins_params+1:]

        assert len(SA) == totcornum, 'error in correlations'
        AR = 0.0  # just to remember to add this
        counter = -1

        for i in range(nins):
            for j in range(len(starflag[i])):
                counter += 1
                passer = -1
                for k in range(ndat):
                    if starflag[i][j] == ins[k]:  #
                        passer += 1
                        COR[i][j][k] = SA[counter] * staract[i][j][passer]

        FMC = 0
        for i in range(len(COR)):
            for j in range(len(COR[i])):
                FMC += COR[i][j]
    else:
        FMC = 0

    MODEL = model(a1, time, kplanets) + offset + ACC + FMC


    for i in range(ndat):
        residuals[i] = rv[i] - MODEL[i]
        for c in range(MOAV):
            if i > c:
                MA = macoef[c][i] * sp.exp(-sp.fabs(time[i-1] - time[i]) / timescale[c][i]) * residuals[i-1]
                residuals[i] -= MA

    inv_sigma2 = 1.0 / (err**2 + jitter**2)
    lnl = sp.sum(residuals ** 2 * inv_sigma2 - sp.log(inv_sigma2)) + sp.log(2*sp.pi) * ndat
    return -0.5 * lnl

def get_at_k1(Ek, Ek_grid, l, cb, conj, phi1k, hphi1k):
    """
    Linearly interpolate the values of phi1k, hphi1k on Ek_grid at point Ek.

    Parameters
    ----------
    Ek : float
        Energy value (not necessarily a grid point).
    Ek_grid : array
        Energy grid.
    l : int
        Lead label.
    cb : int
        Index corresponding to Phi[1](k) matrix element.
    conj : bool
        If conj=True the term in the integral equation is conjugated.
    phi1k : array
        Numpy array with dimensions (len(Ek_grid), nleads, ndm1, ndm0),
        cotaining energy resolved first order density matrix elements Phi[1](k).
    hphi1k : array
        Numpy array with dimensions (len(Ek_grid_ext), nleads, ndm1, ndm0),
        containing Hilbert transform of phi1k.

    Returns
    -------
    array, array
        Interpolated values of phi1k and hphi1k at Ek.
    """
    if Ek<Ek_grid[0] or Ek>Ek_grid[-1]:
        return 0, 0
    #
    #b_idx = np.searchsorted(Ek_grid, Ek)
    b_idx = int((Ek-Ek_grid[0])/(Ek_grid[1]-Ek_grid[0]))+1
    if b_idx == len(Ek_grid): b_idx -= 1
    a_idx = b_idx - 1
    b, a = Ek_grid[b_idx], Ek_grid[a_idx]
    #print(a_idx, b_idx, a, Ek, b)
    #
    fb = phi1k[b_idx, l, cb]
    fa = phi1k[a_idx, l, cb]
    u = (fb-fa)/(b-a)*Ek + (b*fa-a*fb)/(b-a)
    u = u.conjugate() if conj else u
    #
    fb = hphi1k[b_idx, l, cb]
    fa = hphi1k[a_idx, l, cb]
    hu = (fb-fa)/(b-a)*Ek + (b*fa-a*fb)/(b-a)
    hu = hu.conjugate() if conj else hu
    return pi*hu, pi*u