Python scipy.linalg 模块，expm() 实例源码

def entangling_mat(gate):
    """
    Helper function to create the entangling gate matrix
    """
    echoCR = expm(1j * pi / 4 * np.kron(pX, pZ))
    if gate == "CNOT":
        return echoCR
    elif gate == "iSWAP":
        return reduce(lambda x, y: np.dot(y, x),
                      [echoCR, np.kron(C1[6], C1[6]), echoCR])
    elif gate == "SWAP":
        return reduce(lambda x, y: np.dot(y, x),
                      [echoCR, np.kron(C1[6], C1[6]), echoCR, np.kron(
                          np.dot(C1[6], C1[1]), C1[1]), echoCR])
    else:
        raise ValueError("Entangling gate must be one of: CNOT, iSWAP, SWAP.")

def line_search(Y, signs, direction, current_loss, ls_tries):
    '''
    Performs a backtracking line search, starting from Y and W, in the
    direction direction.
    '''
    alpha = 1.
    if current_loss is None:
        current_loss = loss(Y, signs)
    for _ in range(ls_tries):
        Y_new = np.dot(expm(alpha * direction), Y)
        new_loss = loss(Y_new, signs)
        if new_loss < current_loss:
            return True, Y_new, new_loss, alpha
        alpha /= 2.
    else:
        return False, Y_new, new_loss, alpha

def test_expectation():
    """expectation() routine can take a PauliSum operator on a matrix.  Check
    this functionality and the return of the correct scalar value"""
    X = np.array([[0, 1], [1, 0]])

    def RX_gate(phi):
        return expm(-1j*phi*X)

    def rotation_wavefunction(phi):
        state = np.array([[1], [0]])
        return RX_gate(phi).dot(state)

    prog = Program([RX(-2.5)(0)])
    hamiltonian = PauliTerm("Z", 0, 1.0)

    minimizer = MagicMock()
    fake_result = Mock()
    fake_result.fun = 1.0
    minimizer.return_value = fake_result

    fake_qvm = Mock(spec=['wavefunction', 'expectation', 'run'])
    fake_qvm.wavefunction.return_value = (Wavefunction(rotation_wavefunction(-2.5)))
    fake_qvm.expectation.return_value = [0.28366219]
    # for testing expectation
    fake_qvm.run.return_value = [[0], [0]]

    inst = VQE(minimizer)
    energy = inst.expectation(prog, PauliSum([hamiltonian]), None, fake_qvm)
    assert np.isclose(energy, 0.28366219)

    hamiltonian = np.array([[1, 0], [0, -1]])
    energy = inst.expectation(prog, hamiltonian, None, fake_qvm)
    assert np.isclose(energy, 0.28366219)

    prog = Program(H(0))
    hamiltonian = PauliSum([PauliTerm('X', 0)])
    energy = inst.expectation(prog, hamiltonian, 2, fake_qvm)
    assert np.isclose(energy, 1.0)

def rotation_matrix(a):
    """    Creates a 3D rotation matrix for rotation
        around the axis of the vector a. """
    R = eye(4)
    R[:3,:3] = linalg.expm([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])
    return R

def rotation_matrix(a):
    R = eye(4)
    R[:3,:3] = linalg.expm([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])
    return R

def Magnus2(self,direction='x'):
        """Propagate in time using the second order explicit Magnus.
           See: Blanes, Sergio, and Fernando Casas. A concise introduction 
           to geometric numerical integration. Vol. 23. CRC Press, 2016.

           Magnus2 is Eq (4.61), page 128.
        """
        self.reset()
        self.mol.orthoDen()
        self.mol.orthoFock()
        h = -1j*self.stepsize
        for idx,time in enumerate((self.time)):
            if direction.lower() == 'x':
                self.mol.computeDipole()
                self.dipole.append(np.real(self.mol.mu[0]))
            elif direction.lower() == 'y':
                self.mol.computeDipole()
                self.dipole.append(np.real(self.mol.mu[1]))
            elif direction.lower() == 'z':
                self.mol.computeDipole()
                self.dipole.append(np.real(self.mol.mu[2]))

            # record pulse envelope for later plotting, etc.
            self.shape.append(self.pulse(time))
            curDen  = np.copy(self.mol.PO)

            self.addField(time + 0.0*self.stepsize,direction=direction)
            k1 = h*self.mol.FO 
            U = expm(k1)
            self.mol.PO = np.dot(U,np.dot(curDen,self.mol.adj(U))) 
            self.mol.updateFock()

            self.addField(time + 1.0*self.stepsize,direction=direction)
            L  = 0.5*(k1 + h*self.mol.FO)
            U  = expm(L)
            self.mol.PO = np.dot(U,np.dot(curDen,self.mol.adj(U))) 
            self.mol.updateFock()

            # density and Fock are done updating, wrap things up
            self.mol.unOrthoFock()    
            self.mol.unOrthoDen()    
            self.mol.computeEnergy()
            self.Energy.append(np.real(self.mol.energy))

def __init__(self, a_mat, b_vec, settings):
        '''x' = Ax + b  (b_vector here is the constant part of the dynamics, NOT the B input-effect matrix in 'BU') '''

        if isinstance(a_mat, list):
            a_mat = np.array(a_mat, dtype=float)

        if isinstance(b_vec, list):
            b_vec = np.array(b_vec, dtype=float)

        assert isinstance(a_mat, np.ndarray)
        assert isinstance(b_vec, np.ndarray)

        assert isinstance(settings, SimulationSettings)
        assert settings.step > 0

        self.settings = settings
        self.num_dims = b_vec.shape[0]
        assert self.num_dims > 0
        assert a_mat.shape[0] == a_mat.shape[1], "expected A matrix to be a square"
        assert len(b_vec.shape) == 1, "expected passed-in b_vector to be a single row: {}".format(b_vec)
        assert a_mat.shape[0] == b_vec.shape[0], "A matrix and b vector sizes should match"

        self.dy_data = DyData(csr_matrix(a_mat), csr_matrix(b_vec), settings.sparse)

        # initialize simulation result variables
        self.origin_sim = None
        self.vec_values = None
        self.step_offset = None

        # itemsize is bytes per float
        if self.settings.sim_mode == SimulationSettings.SIMULATION:
            mb_per_step = np.dtype(float).itemsize * self.num_dims * self.num_dims / 1024.0 / 1024.0
            self.max_steps_in_mem = max(1, int(settings.sim_in_memory_mb / mb_per_step) - 1)
        elif self.settings.sim_mode == SimulationSettings.MATRIX_EXP:
            self.max_steps_in_mem = 1

            if self.num_dims < 150:
                # dense expm is fast if dims < 150
                self.matrix_exp = dense_expm(a_mat * settings.step).transpose().copy()
            else:
                # sparse expm
                self.matrix_exp = sparse_expm(csc_matrix(a_mat * settings.step)).transpose().toarray()

        self.freeze_attrs()