Python sympy.zeros() Examples

def print_ode(self, latex_output=False):
        '''
        Prints the ode in symbolic form onto the screen/console in actual
        symbols rather than the word of the symbol.

        Parameters
        ----------
        latex_output: bool, optional
            Defaults to false which prints the equation in terms of symbols,
            if set to yes then the formula in terms of latex equations will
            be printed onto the screen.
        '''
        A = self.get_ode_eqn()
        B = sympy.zeros(A.rows,2)
        for i in range(A.shape[0]):
            B[i,0] = sympy.symbols('d' + str(self._stateList[i]) + '/dt=')
            B[i,1] = A[i]

        if latex_output:
            print(sympy.latex(B, mat_str="array", mat_delim=None,
                              inv_trig_style='full'))
        else:
            sympy.pretty_print(B) 

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(3, 1)

        x = sp.Matrix(sp.symbols('x y theta', real=True))
        u = sp.Matrix(sp.symbols('v w', real=True))

        f[0, 0] = u[0, 0] * sp.cos(x[2, 0])
        f[1, 0] = u[0, 0] * sp.sin(x[2, 0])
        f[2, 0] = u[1, 0]

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

def _computeTransitionJacobian(self):
        if self._GMat is None:
            self._computeDependencyMatrix()

        F = sympy.zeros(self.num_transitions, self.num_transitions)
        for i in range(self.num_transitions):
            for j, eqn in enumerate(self._transitionVector):
                for k, state in enumerate(self._iterStateList()):
                    diffEqn = sympy.diff(eqn, state, 1)
                    tempEqn, isDifficult = simplifyEquation(diffEqn)
                    F[i,j] += tempEqn*self._vMat[k,i]
                    self._isDifficult = self._isDifficult or isDifficult

        self._transitionJacobian = F

        self._hasNewTransition.reset('transitionJacobian')
        return F 

def get_equations(self):
        """
        :return: Functions to calculate A, B and f given state x and input u
        """
        f = sp.zeros(6, 1)

        x = sp.Matrix(sp.symbols('rx ry vx vy t w', real=True))
        u = sp.Matrix(sp.symbols('gimbal T', real=True))

        f[0, 0] = x[2, 0]
        f[1, 0] = x[3, 0]
        f[2, 0] = 1 / self.m * sp.sin(x[4, 0] + u[0, 0]) * u[1, 0]
        f[3, 0] = 1 / self.m * (sp.cos(x[4, 0] + u[0, 0]) * u[1, 0] - self.m * self.g)
        f[4, 0] = x[5, 0]
        f[5, 0] = 1 / self.I * (-sp.sin(u[0, 0]) * u[1, 0] * self.r_T)

        f = sp.simplify(f)
        A = sp.simplify(f.jacobian(x))
        B = sp.simplify(f.jacobian(u))

        f_func = sp.lambdify((x, u), f, 'numpy')
        A_func = sp.lambdify((x, u), A, 'numpy')
        B_func = sp.lambdify((x, u), B, 'numpy')

        return f_func, A_func, B_func 

def compute_channel_operation(rho, operators):
        """
        Given a quantum state's density function rho, the effect of the
        channel on this state is:
        rho -> sum_{i=1}^n E_i * rho * E_i^dagger

        Args:
            rho (number): Density function
            operators (list): List of operators

        Returns:
            number: The result of applying the list of operators
        """
        from sympy import zeros
        return sum([E * rho * E.H for E in operators],
                   zeros(operators[0].rows)) 

def get_matrix_from_channel(channel, symbol):
        """
        Extract the numeric parameter matrix.

        Args:
            channel (matrix): a 4x4 symbolic matrix.
            symbol (list): a symbol xi

        Returns:
            matrix: a 4x4 numeric matrix.

        Additional Information:
            Each entry of the 4x4 symbolic input channel matrix is assumed to
            be a polynomial of the form a1x1 + ... + anxn + c. The corresponding
            entry in the output numeric matrix is ai.
        """
        from sympy import Poly
        n = channel.rows
        M = numpy.zeros((n, n), dtype=numpy.complex_)
        for (i, j) in itertools.product(range(n), range(n)):
            M[i, j] = numpy.complex(
                Poly(channel[i, j], symbol).coeff_monomial(symbol))
        return M 

def get_const_matrix_from_channel(channel, symbols):
        """
        Extract the numeric constant matrix.

        Args:
            channel (matrix): a 4x4 symbolic matrix.
            symbols (list): The full list [x1, ..., xn] of symbols
                used in the matrix.

        Returns:
            matrix: a 4x4 numeric matrix.

        Additional Information:
            Each entry of the 4x4 symbolic input channel matrix is assumed to
            be a polynomial of the form a1x1 + ... + anxn + c. The corresponding
            entry in the output numeric matrix is c.
        """
        from sympy import Poly
        n = channel.rows
        M = numpy.zeros((n, n), dtype=numpy.complex_)
        for (i, j) in itertools.product(range(n), range(n)):
            M[i, j] = numpy.complex(
                Poly(channel[i, j], symbols).coeff_monomial(1))
        return M 

def compute_P(self, As):
        """
        This method creates the matrix P in the
        f(x) = 1/2(x*P*x)+q*x
        representation of the objective function
        Args:
            As (list): list of symbolic matrices repersenting the channel matrices

        Returns:
            matrix: The matrix P for the description of the quadaric program
        """
        from sympy import zeros
        vs = [numpy.array(A).flatten() for A in As]
        n = len(vs)
        P = zeros(n, n)
        for (i, j) in itertools.product(range(n), range(n)):
            P[i, j] = 2 * numpy.real(numpy.dot(vs[i], numpy.conj(vs[j])))
        return P 

def _computeTransitionVector(self):
        """
        Get all the transitions into a vector, arranged by state to
        state transition then the birth death processes
        """
        self._transitionVector = sympy.zeros(self.num_transitions, 1)
        _f, _t, eqn = self._unrollTransitionList(self._getAllTransition())
        for i, eqn in enumerate(eqn):
            self._transitionVector[i] = eqn

        return self._transitionVector

    ########################################################################
    #
    # Other type of matrices
    #
    ######################################################################## 

def _computeOdeVector(self):
        # we are only testing it here because we want to be flexible and
        # allow the end user to input more state than initially desired
        if len(self.ode_list) <= self.num_state:
            self._ode = sympy.zeros(self.num_state, 1)
            fromList, _t, eqn = self._unrollTransitionList(self.ode_list)
            for i, eqn in enumerate(eqn):
                if len(fromList[i]) > 1:
                    raise InputError("An explicit ode cannot describe more " +
                                     "than a single state")
                else:
                    self._ode[fromList[i][0]] = eqn
        else:
            raise InputError("The total number of ode is %s " +
                             "where the number of state is %s" %
                             len(self.ode_list), self.num_state)

        return None 

def test_cu3():
    E = eye(2)
    UPPER = Matrix([[1, 0], [0, 0]])
    LOWER = Matrix([[0, 0], [0, 1]])
    theta, phi, lambd = symbols("theta phi lambd")
    U = Circuit().rz(lambd)[0].ry(theta)[0].rz(phi)[0].run_with_sympy_unitary()

    actual_1 = Circuit().cu3(theta, phi, lambd)[0, 1].run(backend="sympy_unitary")
    expected_1 = reduce(TensorProduct, [E, UPPER]) + reduce(TensorProduct, [U, LOWER])
    print("actual")
    print(simplify(actual_1))
    print("expected")
    print(simplify(expected_1))
    print("diff")
    print(simplify(actual_1 - expected_1))
    assert simplify(actual_1 - expected_1) == zeros(4) 

def test_rygate():
    t = symbols('t')
    u = Circuit().ry(t)[0].to_unitary()
    expected = simplify(exp(-I * t / 2 * Circuit().y[0].to_unitary()))
    assert simplify(u - expected) == zeros(2) 

def FdiagPlus1(a,n):
    f = Fdiag(a,n-1)
    f = f.col_insert(n-1, zeros(n-1,1))
    f = f.row_insert(n-1, Matrix(1,n, lambda i,j: (1 if j==n-1 else 0)))
    return f 

def test_u1():
    lambd = symbols("lambd")

    actual_1 = Circuit().u1(lambd)[0].run(backend="sympy_unitary")
    expected_1 = Circuit().rz(lambd)[0].run_with_sympy_unitary()
    assert simplify(actual_1 - expected_1) == zeros(2) 

def test_u2():
    phi, lambd = symbols("phi lambd")

    actual_1 = Circuit().u2(phi, lambd)[0].run(backend="sympy_unitary")
    expected_1 = Circuit().rz(lambd)[0].ry(pi / 2)[0].rz(phi)[0].run_with_sympy_unitary()
    assert simplify(actual_1 - expected_1) == zeros(2) 

def test_u3():
    theta, phi, lambd = symbols("theta phi lambd")

    actual_1 = Circuit().u3(theta, phi, lambd)[0].run(backend="sympy_unitary")
    assert actual_1[0, 0] != 0
    expected_1 = Circuit().rz(lambd)[0].ry(theta)[0].rz(phi)[0].run_with_sympy_unitary()
    assert expected_1[0, 0] != 0
    assert simplify(actual_1 - expected_1) == zeros(2) 

def test_cu1():
    E = eye(2)
    UPPER = Matrix([[1, 0], [0, 0]])
    LOWER = Matrix([[0, 0], [0, 1]])
    lambd = symbols("lambd")
    U = Circuit().rz(lambd)[0].run_with_sympy_unitary()
    U /= U[0, 0]

    actual_1 = Circuit().cu1(lambd)[0, 1].run(backend="sympy_unitary")
    actual_1 /= actual_1[0, 0]
    expected_1 = reduce(TensorProduct, [UPPER, E]) + reduce(TensorProduct, [LOWER, U])
    assert simplify(actual_1 - expected_1) == zeros(4) 

def test_cygate():
    u1 = Circuit().cy[1, 0].to_unitary()
    u2 = Matrix([
        [1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 0, -I],
        [0, 0, I, 0]])
    assert simplify(u1 - u2) == zeros(4) 

def test_rxgate():
    t = symbols('t')
    u = Circuit().rx(t)[0].to_unitary()
    expected = exp(-I * t / 2 * Circuit().x[0].to_unitary())
    assert simplify(u - expected) == zeros(2) 

def _computeDependencyMatrix(self):
        """
        Obtain the dependency matrix/graph. G_{i,j} indicate whether invoking
        the transition j will cause the rate to change for transition j
        """
        if self._lambdaMat is None:
            self._computeReactantMatrix()
        if self._vMat is None:
            self._computeStateChangeMatrix()

        nt = self.num_transitions
        self._GMat = np.zeros((nt, nt), int)

        for i in range(nt):
            for j in range(nt):
                d = 0
                for k in range(self.num_state):
                    d = d or (self._lambdaMat[k, i] and self._vMat[k, j])
                self._GMat[i, j] = d

        return self._GMat


    ########################################################################
    # Unrolling of the information
    # state
    ######################################################################## 

def test_rzgate():
    t = symbols('t')
    u = Circuit().rz(t)[0].to_unitary()
    expected = exp(-I * t / 2 * Circuit().z[0].to_unitary())
    assert simplify(u - expected) == zeros(2) 

def test_rxxgate():
    t = symbols('t')
    u = Circuit().rxx(t)[1, 0].to_unitary()
    expected = simplify(exp(-I * t / 2 * Circuit().x[0, 1].to_unitary()))
    assert simplify(u - expected) == zeros(4) 

def test_rzzgate():
    t = symbols('t')
    u = Circuit().rzz(t)[1, 0].to_unitary()
    expected = simplify(exp(-I * t / 2 * Circuit().z[0, 1].to_unitary()))
    assert simplify(u - expected) == zeros(4) 

def test_swapgate():
    assert simplify(Circuit().swap[1, 0].to_unitary() -
                    Matrix([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])) == zeros(4) 

def test_cswapgate():
    expected = eye(8)
    expected[4:, 4:] = Circuit().swap[1, 0].to_unitary()
    assert simplify(Circuit().cswap[2, 1, 0].to_unitary() - expected) == zeros(8) 

def _computeBirthDeathRate(self):
        '''
        Note that this is different to _birthDeathVector because
        this is of length (number of birth and death process) while
        _birthDeathVector in baseOdeModel has the same length as
        the number of states
        '''
        if self.num_birth_deaths == 0:
            A = sympy.zeros(1, 1)
        else:
            A = sympy.zeros(self.num_birth_deaths, 1)

            # go through all the transition objects
            for i, bd in enumerate(self.birth_death_list):
                A[i] += eval(self._checkEquation(bd.equation()))

        # assign back
        self._birthDeathRate = A
        # compilation of the symbolic calculation.  Note here that we are
        # going to recompile the total transitions because it might
        # have changed
        f = self._SC.compileExprAndFormat
        if self._isDifficult:
            self._birthDeathRateCompile = f(self._sp,
                                            self._birthDeathRate,
                                            modules='mpmath')
        else:
            self._birthDeathRateCompile = f(self._sp,
                                            self._birthDeathRate)

        self._hasNewTransition.reset('birthDeathRateCompile')

        return None 

def generateDirectedDependencyGraph(ode_matrix, transition=None):
    """
    Returns a binary matrix that contains the direction of the transition in
    a state

    Parameters
    ----------
    ode_matrix: :class:`sympy.matrcies.MatrixBase`
        A matrix of size [number of states x 1].  Obtained by
        invoking :meth:`DeterministicOde.get_ode_eqn`
    transition: list, optional
        list of transitions.  Can be generated by
        :func:`getMatchingExpressionVector`

    Returns
    -------
    G: :class:`numpy.ndarray`
        Two dimensional array of size [number of state x number of transitions]
        where each column has two entry,
        -1 and 1 to indicate the direction of the transition and the state.
        All column sum to one, i.e. transition must have a source and target.
    """
    assert isinstance(ode_matrix, MatrixBase), \
        "Expecting a vector of expressions"

    if transition is None:
        transition = getMatchingExpressionVector(ode_matrix, True)
    else:
        assert isinstance(transition, list), "Require a list of transitions"

    B = np.zeros((len(ode_matrix), len(transition)))
    for i, a in enumerate(ode_matrix):
        for j, transitionTuple in enumerate(transition):
            t1, t2 = transitionTuple
            if _hasExpression(a, t1):
                B[i, j] += -1  # going out
            if _hasExpression(a, t2):
                B[i, j] += 1   # coming in
    return B 

def _extractUpperTriangle(self, A, nrow=None, ncol=None):
        """
        Extract the upper triangle of matrix A

        Parameters
        ----------
        A: :mod:`sympy.matrices.matrices`
            input matrix
        nrow: int
            number of row
        ncol: int
            number of column

        Returns
        -------
        :mod:`sympy.matrices.matrices`
            An upper triangle matrix

        """
        if nrow is None:
            nrow = len(A[:, 0])

        if ncol is None:
            ncol = len(A[0, :])

        B = sympy.zeros(nrow, ncol)
        for i in range(0, nrow):
            for j in range(i, ncol):
                B[i,j] = A[i, j]

        return B 

def _computeStateChangeMatrix(self):
        """
        The state change matrix, where
        .. math::
            v_{i,j} = \\left\\{ 1, &if transition j cause state i to lose a particle, \\\\
                               -1, &if transition j cause state i to gain a particle, \\\\
                                0, &otherwise \\right.
        """
        self._vMat = np.zeros((self.num_state, self.num_transitions), int)

        f, t, eqn = self._unrollTransitionList(self._getAllTransition())
        for j, _eqn in enumerate(eqn):
            if j < self.num_pure_transitions:
                for k1 in f[j]:
                    self._vMat[k1, j] += -1
                for k2 in t[j]:
                    self._vMat[k2, j] += 1
            else:
                bdObj = self._birthDeathList[j - self.num_pure_transitions]
                if bdObj.transition_type is TransitionType.B:
                    for k1 in f[j]:
                        self._vMat[k1, j] += 1
                elif bdObj.transition_type is TransitionType.D:
                    for k2 in f[j]:
                        self._vMat[k2, j] += -1

        return self._vMat 

def _computeBirthDeathVector(self):
        # holder
        self._birthDeathVector = sympy.zeros(self.num_state, 1)
        # go through all the transition objects
        for bdObj in self._birthDeathList:
            fromIndex, _toIndex, eqn = self._unrollTransition(bdObj)
            for i in fromIndex:
                if bdObj.transition_type is TransitionType.B:
                    self._birthDeathVector[i] += eqn
                elif bdObj.transition_type is TransitionType.D:
                    self._birthDeathVector[i] -= eqn

        return self._birthDeathVector