Python linalg.expm函数代码示例

def align_magnetism(m, vectors):
    """ Rotates a matrix, to align its components with the direction
  of the magnetism """
    if not len(m) == 2 * len(vectors):  # stop if they don't have
        # compatible dimensions
        raise
    # pauli matrices
    from scipy.sparse import csc_matrix, bmat

    sx = csc_matrix([[0.0, 1.0], [1.0, 0.0]])
    sy = csc_matrix([[0.0, -1j], [1j, 0.0]])
    sz = csc_matrix([[1.0, 0.0], [0.0, -1.0]])
    n = len(m) / 2  # number of sites
    R = [[None for i in range(n)] for j in range(n)]  # rotation matrix
    from scipy.linalg import expm  # exponenciate matrix

    for (i, v) in zip(range(n), vectors):  # loop over sites
        vv = np.sqrt(v.dot(v))  # norm of v
        if vv > 0.000001:  # if nonzero scale
            u = v / vv
        else:  # if zero put to zero
            u = np.array([0.0, 0.0, 0.0])
        #    rot = u[0]*sx + u[1]*sy + u[2]*sz
        uxy = np.sqrt(u[0] ** 2 + u[1] ** 2)  # component in xy plane
        phi = np.arctan2(u[1], u[0])
        theta = np.arctan2(uxy, u[2])
        r1 = phi * sz / 2.0  # rotate along z
        r2 = theta * sy / 2.0  # rotate along y
        # a factor 2 is taken out due to 1/2 of S
        rot = expm(1j * r2) * expm(1j * r1)
        R[i][i] = rot  # save term
    R = bmat(R)  # convert to full sparse matrix
    mout = R * csc_matrix(m) * R.H  # rotate matrix
    return mout.todense()  # return dense matrix

def make_random_su(N):
    
    A = spla.expm(2 * np.pi * 1j * npr.random((N, N)))
    B = (A - np.trace(A) * np.identity(N) / N)
    C = 0.5 * (B - np.conj(B.T))
    
    return spla.expm(C)

    def _solve_numericaly(self, u, x0, t, t_list, t0, dps=2):
        """ returns the numeric evaluation of the system for input u, know state x0 at time t0 and times t_list
        """
        result = []
        for t_i in t_list:
            # we use the arbitrary precision module mpmath for numercial evaluation of the matrix exponentials
            first = (
                np.array(np.array(self.represent[2]), np.float)
                .dot(expm(np.array(np.array((self.represent[0] * (t_i - t0)).evalf()), np.float)))
                .dot(np.array(np.array(x0), np.float))
            )

            second = np.array(np.array((self.represent[3] * u.subs(t, t_i)).evalf()), np.float)

            integrand = (
                lambda tau: np.array(np.array(self.represent[2]), np.float)
                .dot(expm(np.array(np.array((self.represent[0] * (t_i - tau)).evalf()), np.float)))
                .dot(np.array(np.array(self.represent[1]), np.float))
                .dot(np.array(np.array(u.subs(t, tau).evalf()), np.float))
            )

            # the result must have the same shape as D:
            integral = zeros(self.represent[2].rows, 1)

            # Loop through every entry and evaluate the integral using mpmath.quad()
            for row_idx in xrange(self.represent[2].rows):

                integral[row_idx, 0] = quad(lambda x: integrand(x)[row_idx, 0], t0, t_i)[0]

            result.append(Matrix(first) + Matrix(second) + integral)

        # return sum of results
        return result

    def make_propagators(pb=0.0, kex=0.0, dw=0.0, r_coxy=5.0, dr_coxy=0.0,
                         r_nz=1.5, r_2coznz=0.0, etaxy=0.0, etaz=0.0,
                         j_nco=0.0, dj_nco=0.0, cs_offset=0.0):

        w1 = 2.0 * pi / (4.0 * pwco90)
        l_free, l_w1x, l_w1y = compute_liouvillians(pb=pb, kex=kex, dw=dw,
                                                    r_coxy=r_coxy, dr_coxy=dr_coxy,
                                                    r_nz=r_nz, r_2coznz=r_2coznz,
                                                    etaxy=etaxy, etaz=etaz,
                                                    j_nco=j_nco, dj_nco=dj_nco,
                                                    cs_offset=cs_offset, w1=w1)

        p_equil = expm(l_free * time_equil)
        p_taucc = expm(l_free * taucc)
        p_neg = expm(l_free * -2.0 * pwco90 / pi)
        p_90py = expm((l_free + l_w1y) * pwco90)
        p_90my = expm((l_free - l_w1y) * pwco90)
        p_180px = P180X  # Perfect 180 for CPMG blocks
        p_180py = matrix_power(p_90py, 2)
        p_180my = matrix_power(p_90py, 2)

        ps = (p_equil, p_taucc, p_neg, p_90py, p_90my,
              p_180px, p_180py, p_180my)

        return l_free, ps

    def _learnStep(self):
        """ Main part of the algorithm. """
        I = eye(self.numParameters)
        self._produceSamples()
        utilities = self.shapingFunction(self._currentEvaluations)
        utilities /= sum(utilities)  # make the utilities sum to 1
        if self.uniformBaseline:
            utilities -= 1./self.batchSize
        samples = array(map(self._base2sample, self._population))

        dCenter = dot(samples.T, utilities)
        covGradient = dot(array([outer(s,s) - I for s in samples]).T, utilities)
        covTrace = trace(covGradient)
        covGradient -= covTrace/self.numParameters * I
        dA = 0.5 * (self.scaleLearningRate * covTrace/self.numParameters * I
                    +self.covLearningRate * covGradient)

        self._lastLogDetA = self._logDetA
        self._lastInvA = self._invA

        self._center += self.centerLearningRate * dot(self._A, dCenter)
        self._A = dot(self._A, expm(dA))
        self._invA = dot(expm(-dA), self._invA)
        self._logDetA += 0.5 * self.scaleLearningRate * covTrace
        if self.storeAllDistributions:
            self._allDistributions.append((self._center.copy(), self._A.copy()))

def transfer_f(dw,aas,aai,eps,deltaw,f):
    """
    Args:
    dw: size of the grid spacing
    aas=relative slowness of the signal mode
    aai=relative slowness of the idler mode
    lnl=inverse of the strength of the nonlinearity
    deltaw:  specifies the size of the frequency grid going from
    -deltaw to deltaw for each frequency
    f: shape of the pump function
    """
    ddws=np.arange(-deltaw-dw/2,deltaw+dw/2,dw)
    deltaks=aas*ddws
    ddwi=np.arange(-deltaw-dw/2,deltaw+dw/2,dw)
    deltaki=aai*ddwi
    ds=np.diag(deltaks)
    di=np.diag(deltaki)


    def ff(x,y):
        return f(x+y)
    
    v=eps*(dw)*ff(ddwi[:,None],ddws[None,:])
    G=1j*np.concatenate((np.concatenate((ds,v),axis=1),np.concatenate((-v,-di),axis=1)),axis=0)
    z=1;
    dsi=np.concatenate((deltaks,-deltaki),axis=0)
    U0=linalg.expm(-1j*np.diag(dsi)*z/2)
    GG=np.dot(np.dot(U0,linalg.expm(G)),U0)
    n=len(ddws)
    return (GG[0:n,0:n],GG[n:2*n,0:n],GG[0:n,n:2*n],GG[n:2*n,n:2*n])

def ghz_simult_trajectory(stages):
    rhos = []
    for stage in stages:
        U = np.eye(8, dtype=complex)
        I = pyle.tomo.sigmaI
        X = pyle.tomo.sigmaX
        Y = pyle.tomo.sigmaY
        couple = (pyle.tensor((X,X,I)) + pyle.tensor((Y,Y,I)) +
                  pyle.tensor((X,I,X)) + pyle.tensor((Y,I,Y)) +
                  pyle.tensor((I,X,X)) + pyle.tensor((I,Y,Y))) / 2.0
        
        if stage > 0:
            fraction = np.clip(stage-0, 0, 1)
            H = -1j * (np.pi/2)/2 * (pyle.tensor((Y,I,I)) + pyle.tensor((I,Y,I)) + pyle.tensor((I,I,Y)))
            U = np.dot(expm(fraction * H), U)
    
        if stage > 1:
            fraction = np.clip(stage-1, 0, 1)
            H = -1j * np.pi/2 * couple
            U = np.dot(expm(fraction * H), U)
        
        if stage > 2:
            fraction = np.clip(stage-2, 0, 1)
            H = -1j * (np.pi/2)/2 * (pyle.tensor((X,I,I)) + pyle.tensor((I,X,I)) + pyle.tensor((I,I,X)))
            U = np.dot(expm(fraction * H), U)
    
        psi0 = np.array([1,0,0,0,0,0,0,0], dtype=complex)
        psi_th = np.dot(U, psi0)
        rho_th = pyle.ket2rho(psi_th)
        rhos.append(rho_th)
    return np.array(rhos)

def linear_ode_discretation(F, L=None, Q=None, dt=1):
    n = F.shape[0]

    if L is None:
        L = eye(n)

    if Q is None:
        Q = zeros((n,n))

    A = expm(F*dt)

    phi = zeros((2*n, 2*n))

    phi[0:n,     0:n] = F
    phi[0:n,   n:2*n] = L.dot(Q).dot(L.T)
    phi[n:2*n, n:2*n] = -F.T

    zo = vstack((zeros((n,n)), eye(n)))

    CD = expm(phi*dt).dot(zo)

    C = CD[0:n,:]
    D = CD[n:2*n,:]
    q = C.dot(inv(D))

    return (A, q)

def Iterate_LS_Strain(xys0_pix, xys1_pix, num = 10):
    solve = Get_LS_DGradient(xys0_pix, xys1_pix)
    for j in range(0,num):
        xy_pix = Simulate_Shifted_Spots(xys0_pix, linalg.expm(solve), det2lab_mat)
        dsolve = Get_LS_DGradient(xy_pix, xys1_pix)
        solve = solve+dsolve
    return linalg.expm(solve)

def Char_Gate(NV,res ,B_field=400):
    """
    Characterize the gate, take the NV centre, the resonance paramters and the Bfield as input
    returns the fidelity with which an x-gate can be implemented.
    """


    #data = np.loadtxt("NV_Sim_8.dat") #Placeholder data to test the script
    #NV = np.vstack((data[:,3],data[:,4]))
    #physical constants
    gamma_c = 1.071e3 #g-factor for C13 in Hz/G
    #Model parameters
    omega_larmor = 2*np.pi*gamma_c*B_field
    tau_larmor = 2*np.pi/omega_larmor
    tau = res[0]
    n_pulses = int(res[1]*2) #So that we do a pi -pulse

    Ix = 0.5 * np.array([[0,1],[1,0]])
    Iz = 0.5* np.array([[1,0],[0,-1]])
    H0 = (omega_larmor)*Iz
    exH0 =linalg.expm(-1j*H0*tau)

    S_final =1
    for idC in range(np.shape(NV)[1]):
        A= 2*np.pi*NV[0,idC]
        B= 2*np.pi*NV[1,idC]  #Converts to radial frequency in Hz/G
        H1 = (A+omega_larmor) *Iz +B*Ix
        exH1 = linalg.expm(-1j*H1*tau)
        V0 = exH0.dot(exH1.dot(exH1.dot(exH0)))
        V1 = exH1.dot(exH0.dot(exH0.dot(exH1)))
        S = np.real(np.trace(np.dot(np.linalg.matrix_power(V0,n_pulses),np.linalg.matrix_power(V1,n_pulses)))/2)
        S_final = S_final *S
    F = (1-(S_final+1)/2) #Converting from probability of measuring +X to fidelity of -X (x-rotation)
    return F

def discretize(F,G,Q,Ts):
    Phi = sp_linalg.expm(F*Ts)

    # Control matrix Lambda, not to be used
    L = np.zeros((F.shape[0],1))

    A_z = np.zeros((L.shape[1], F.shape[1]+L.shape[1]))
    A1 = np.vstack((np.hstack((F,L)),A_z))

    Loan1 = sp_linalg.expm(A1*Ts)
    Lambda = Loan1[0:L.shape[0], F.shape[1]:F.shape[1]+L.shape[1]]

    # Covariance
    Qc = symmetrize(np.dot(G, np.dot(Q, G.T)))
    dim = F.shape[0]
    A2 = np.vstack((np.hstack((-F, Qc)), np.hstack((np.zeros((dim,dim)), F.T))))
    Loan2 = sp_linalg.expm(A2*Ts)
    G2 = Loan2[0:dim, dim:2*dim]
    F3 = Loan2[dim:2*dim, dim:2*dim]

    # Calculate Gamma*Gamma.T
    Qd = symmetrize(np.dot(F3.T, G2))
    L = np.linalg.cholesky(Qd)
    
    return Phi, Qd, L

    def pdf(self, x):
        """ probability density function """
        if not np.isscalar(x):
            x = np.asarray(x)
            res = np.zeros_like(x)
            nz = (x > 0)
            if np.any(nz):
                if self.method == 'sum':
                    factor = np.exp(-x[nz, None] * self.rates[..., :]) \
                                / self.rates[..., :]
                    res[nz] = np.sum(self._terms[..., :] * factor, axis=1)
                else:
                    Theta = (np.diag(-self.rates, 0) + 
                             np.diag(self.rates[:-1], 1))
                    for i in np.flatnonzero(nz):
                        res.flat[i] = \
                                1 - linalg.expm(x.flat[i]*Theta)[0, :].sum()
 
        elif x == 0:
            res = 0
        else:
            if self.method == 'sum':
                factor = np.exp(-x*self.rates)/self.ratesx
                res[nz] = np.sum(self._terms * factor)
            else:
                Theta = np.diag(-self.rates, 0) + np.diag(self.rates[:-1], 1)
                res = 1 - linalg.expm(x*Theta)[0, :].sum()
        return res

 def test_padecases_float(self):
     # test single-precision cases
     a1 = eye(3, dtype=float)*1e-1; e1 = exp(1e-1)*eye(3)
     a2 = eye(3, dtype=float);      e2 = exp(1.0)*eye(3)
     a3 = eye(3, dtype=float)*10;   e3 = exp(10.)*eye(3)
     assert_array_almost_equal(expm(a1),e1)
     assert_array_almost_equal(expm(a2),e2)
     assert_array_almost_equal(expm(a3),e3)

    def test_consistency(self):
        a = array([[0.,1],[-1,0]])
        assert_array_almost_equal(expm(a), expm2(a))
        assert_array_almost_equal(expm(a), expm3(a))

        a = array([[1j,1],[-1,-2j]])
        assert_array_almost_equal(expm(a), expm2(a))
        assert_array_almost_equal(expm(a), expm3(a))

def test_QobjExpmExplicitDense():
    "Qobj expm (explicit dense)"
    data = np.random.random(
        (15, 15)) + 1j * np.random.random((15, 15)) - (0.5 + 0.5j)
    A = Qobj(data)
    B = A.expm(method='dense')
    assert_((B.data.todense() - np.matrix(la.expm(data)) < 1e-10).all())
    B = A.expm(method='scipy-delse')
    assert_((B.data.todense() - np.matrix(la.expm(data)) < 1e-10).all())