我们从Python开源项目中，提取了以下22个代码示例，用于说明如何使用sympy.zeros()。
def test_Dirac(): gamma0 = mgamma(0) gamma1 = mgamma(1) gamma2 = mgamma(2) gamma3 = mgamma(3) gamma5 = mgamma(5) # gamma*I -> I*gamma (see #354) assert gamma5 == gamma0 * gamma1 * gamma2 * gamma3 * I assert gamma1 * gamma2 + gamma2 * gamma1 == zeros(4) assert gamma0 * gamma0 == eye(4) * minkowski_tensor[0, 0] assert gamma2 * gamma2 != eye(4) * minkowski_tensor[0, 0] assert gamma2 * gamma2 == eye(4) * minkowski_tensor[2, 2] assert mgamma(5, True) == \ mgamma(0, True)*mgamma(1, True)*mgamma(2, True)*mgamma(3, True)*I
def mass_matrix_full(self): """ Augments the coefficients of qdots to the mass_matrix. """ n = len(self._q) if self.eom is None: raise ValueError('Need to compute the equations of motion first') #THE FIRST TWO ROWS OF THE MATRIX row1 = eye(n).row_join(zeros(n, n)) row2 = zeros(n, n).row_join(self.mass_matrix) if self.coneqs is not None: m = len(self.coneqs) I = eye(n).row_join(zeros(n, n + m)) below_eye = zeros(n + m, n) A = (self.mass_matrix).col_join((self._m_cd).row_join(zeros(m, m))) below_I = below_eye.row_join(A) return I.col_join(below_I) else: A = row1.col_join(row2) return A
def dsr_graph_adj(self, keep = None): """Return the adjacency matrix of the directed species-reaction graph. Optionally remove the variables of the influence matrix not in *keep*. References: Feliu, E., & Wiuf, C. (2015). Finding the positive feedback loops underlying multi-stationarity. BMC systems biology, 9(1), 1 """ A = self.stoich_matrix im = self.influence_matrix() if keep is not None: im = im.applyfunc(lambda x: x if x in keep else 0) im = im.applyfunc(lambda x: -1 if str(x)[0] == "-" else (1 if x != 0 else 0)) adjacency = sp.zeros(im.rows, im.rows).row_join(im).col_join(A.T.row_join(sp.zeros(A.cols, A.cols))) return adjacency ### Connectivity properties ###
def getRawVelocities(self, subs={}, include_fast_reaction=True, include_slow_reaction=True): definition = self.__definition.getDeveloppedInternalMathFormula() if not KineticLawIdentifier.isReversible(self, definition): res = zeros(1) if (self.reaction.fast and include_fast_reaction) or (not self.reaction.fast and include_slow_reaction): res[0] += definition.subs(subs) return res else: res_front = zeros(1) res_backward = zeros(1) if (self.reaction.fast and include_fast_reaction) or (not self.reaction.fast and include_slow_reaction): (definition_front, definition_back) = KineticLawIdentifier.getReversibleRates(self, definition) res_front[0] += definition_front.subs(subs) res_backward[0] += definition_back.subs(subs) return res_front.col_join(res_backward)
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 _compute_basis(self, closed_form): """ Compute a set of functions B=(f1, ..., fn), a nxn matrix M and a 1xn matrix C such that: closed_form = C B z d/dz B = M B. """ from sympy.matrices import Matrix, eye, zeros afactors = [_x + a for a in self.func.ap] bfactors = [_x + b - 1 for b in self.func.bq] expr = _x*Mul(*bfactors) - self.z*Mul(*afactors) poly = Poly(expr, _x) n = poly.degree() - 1 b = [closed_form] for _ in xrange(n): b.append(self.z*b[-1].diff(self.z)) self.B = Matrix(b) self.C = Matrix([[1] + [0]*n]) m = eye(n) m = m.col_insert(0, zeros(n, 1)) l = poly.all_coeffs()[1:] l.reverse() self.M = m.row_insert(n, -Matrix([l])/poly.all_coeffs()[0])
def test_dcm(): q1, q2, q3, q4 = dynamicsymbols('q1 q2 q3 q4') N = ReferenceFrame('N') A = N.orientnew('A', 'Axis', [q1, N.z]) B = A.orientnew('B', 'Axis', [q2, A.x]) C = B.orientnew('C', 'Axis', [q3, B.y]) D = N.orientnew('D', 'Axis', [q4, N.y]) E = N.orientnew('E', 'Space', [q1, q2, q3], '123') assert N.dcm(C) == Matrix([ [- sin(q1) * sin(q2) * sin(q3) + cos(q1) * cos(q3), - sin(q1) * cos(q2), sin(q1) * sin(q2) * cos(q3) + sin(q3) * cos(q1)], [sin(q1) * cos(q3) + sin(q2) * sin(q3) * cos(q1), cos(q1) * cos(q2), sin(q1) * sin(q3) - sin(q2) * cos(q1) * cos(q3)], [- sin(q3) * cos(q2), sin(q2), cos(q2) * cos(q3)]]) # This is a little touchy. Is it ok to use simplify in assert? test_mat = D.dcm(C) - Matrix( [[cos(q1) * cos(q3) * cos(q4) - sin(q3) * (- sin(q4) * cos(q2) + sin(q1) * sin(q2) * cos(q4)), - sin(q2) * sin(q4) - sin(q1) * cos(q2) * cos(q4), sin(q3) * cos(q1) * cos(q4) + cos(q3) * (- sin(q4) * cos(q2) + sin(q1) * sin(q2) * cos(q4))], [sin(q1) * cos(q3) + sin(q2) * sin(q3) * cos(q1), cos(q1) * cos(q2), sin(q1) * sin(q3) - sin(q2) * cos(q1) * cos(q3)], [sin(q4) * cos(q1) * cos(q3) - sin(q3) * (cos(q2) * cos(q4) + sin(q1) * sin(q2) * sin(q4)), sin(q2) * cos(q4) - sin(q1) * sin(q4) * cos(q2), sin(q3) * sin(q4) * cos(q1) + cos(q3) * (cos(q2) * cos(q4) + sin(q1) * sin(q2) * sin(q4))]]) assert test_mat.expand() == zeros(3, 3) assert E.dcm(N) == Matrix( [[cos(q2)*cos(q3), sin(q3)*cos(q2), -sin(q2)], [sin(q1)*sin(q2)*cos(q3) - sin(q3)*cos(q1), sin(q1)*sin(q2)*sin(q3) + cos(q1)*cos(q3), sin(q1)*cos(q2)], [sin(q1)*sin(q3) + sin(q2)*cos(q1)*cos(q3), - sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1), cos(q1)*cos(q2)]])
def test_function_return_types(): # Lets ensure that decompositions of immutable matrices remain immutable # I.e. do MatrixBase methods return the correct class? X = ImmutableMatrix([[1, 2], [3, 4]]) Y = ImmutableMatrix([[1], [0]]) q, r = X.QRdecomposition() assert (type(q), type(r)) == (ImmutableMatrix, ImmutableMatrix) assert type(X.LUsolve(Y)) == ImmutableMatrix assert type(X.QRsolve(Y)) == ImmutableMatrix X = ImmutableMatrix([[1, 2], [2, 1]]) assert X.T == X assert X.is_symmetric assert type(X.cholesky()) == ImmutableMatrix L, D = X.LDLdecomposition() assert (type(L), type(D)) == (ImmutableMatrix, ImmutableMatrix) assert X.is_diagonalizable() assert X.berkowitz_det() == -3 assert X.norm(2) == 3 assert type(X.eigenvects()[0][2][0]) == ImmutableMatrix assert type(zeros(3, 3).as_immutable().nullspace()[0]) == ImmutableMatrix X = ImmutableMatrix([[1, 0], [2, 1]]) assert type(X.lower_triangular_solve(Y)) == ImmutableMatrix assert type(X.T.upper_triangular_solve(Y)) == ImmutableMatrix assert type(X.minorMatrix(0, 0)) == ImmutableMatrix # issue 6279 # https://github.com/sympy/sympy/issues/6279 # Test that Immutable _op_ Immutable => Immutable and not MatExpr
def vandermonde(order, dim=1, syms='a b c d'): """Comptutes a Vandermonde matrix of given order and dimension. Define syms to give beginning strings for temporary variables. Returns the Matrix, the temporary variables, and the terms for the polynomials. """ syms = syms.split() n = len(syms) if n < dim: new_syms = [] for i in range(dim - n): j, rem = divmod(i, n) new_syms.append(syms[rem] + str(j)) syms.extend(new_syms) terms = [] for i in range(order + 1): terms.extend(comb_w_rep(dim, i)) rank = len(terms) V = zeros(rank) generators = [symbol_gen(syms[i]) for i in range(dim)] all_syms = [] for i in range(rank): row_syms = [next(g) for g in generators] all_syms.append(row_syms) for j, term in enumerate(terms): v_entry = 1 for k in term: v_entry *= row_syms[k] V[i*rank + j] = v_entry return V, all_syms, terms
def zassenhaus(basis_a, basis_b): r""" Given two bases ``basis_a`` and ``basis_b`` of vector spaces :math:`U, W \subseteq V`, computes bases of :math:`U + W` and :math:`U \cap W` using the Zassenhaus algorithm. """ # handle the case where one of the bases is empty if len(basis_a) == 0 or len(basis_b) == 0: mat, pivot = sp.Matrix(basis_a).rref() plus_basis = mat[:len(pivot), :].tolist() return ZassenhausResult(sum=plus_basis, intersection=[]) else: # Set up the Zassenhaus matrix from the given bases. A = sp.Matrix(basis_a) B = sp.Matrix(basis_b) dim = A.shape[1] if B.shape[1] != dim: raise ValueError('Inconsistent dimensions of the two bases given.') zassenhaus_mat = A.row_join(A).col_join(B.row_join(sp.zeros(*B.shape))) mat, pivot = zassenhaus_mat.rref() # idx is the row index of the first row belonging to the intersection basis idx = np.searchsorted(pivot, dim) plus_basis = mat[:idx, :dim].tolist() # The length of the pivot table is used to get rid of all-zero rows int_basis = mat[idx:len(pivot), dim:].tolist() return ZassenhausResult(sum=plus_basis, intersection=int_basis)
def hermitian_to_vector(matrix, basis): """ Returns a the vector representing the ``matrix`` w.r.t. the given *orthonormal* ``basis``. """ _assert_orthogonal(basis) vec = tuple( frobenius_product(matrix, b) / frobenius_product(b, b) for b in basis ) vec = tuple(v.nsimplify() for v in vec) # check consistency assert matrix.equals( sum((v * b for v, b in zip(vec, basis)), sp.zeros(*matrix.shape)) ) return vec
def _compute_basis(self, closed_form): """ Compute a set of functions B=(f1, ..., fn), a nxn matrix M and a 1xn matrix C such that: closed_form = C B z d/dz B = M B. """ from sympy.matrices import Matrix, eye, zeros afactors = [_x + a for a in self.func.ap] bfactors = [_x + b - 1 for b in self.func.bq] expr = _x*Mul(*bfactors) - self.z*Mul(*afactors) poly = Poly(expr, _x) n = poly.degree() - 1 b = [closed_form] for _ in range(n): b.append(self.z*b[-1].diff(self.z)) self.B = Matrix(b) self.C = Matrix([[1] + [0]*n]) m = eye(n) m = m.col_insert(0, zeros(n, 1)) l = poly.all_coeffs()[1:] l.reverse() self.M = m.row_insert(n, -Matrix([l])/poly.all_coeffs()[0])
def is_ss_flux(self, v): """Check if v is a steady state flux.""" if len(v) != self.n_reactions: raise ValueError("Invalid flux vector: wrong length.") return self.stoich_matrix * sp.Matrix(v) == sp.zeros(self.n_species, 1)
def is_cyclic_ss_flux(self, v): """Check if v is a cyclic steady state flux.""" if len(v) != self.n_reactions: raise ValueError("Invalid flux vector: wrong length.") return self.incidence_matrix * sp.Matrix(v) == sp.zeros(self.n_complexes, 1)
def build_hypergeometric_formula(func): """ Create a formula object representing the hypergeometric function ``func``. """ # We know that no `ap` are negative integers, otherwise "detect poly" # would have kicked in. However, `ap` could be empty. In this case we can # use a different basis. # I'm not aware of a basis that works in all cases. from sympy import zeros, Matrix, eye z = Dummy('z') if func.ap: afactors = [_x + a for a in func.ap] bfactors = [_x + b - 1 for b in func.bq] expr = _x*Mul(*bfactors) - z*Mul(*afactors) poly = Poly(expr, _x) n = poly.degree() basis = [] M = zeros(n) for k in xrange(n): a = func.ap[0] + k basis += [hyper([a] + list(func.ap[1:]), func.bq, z)] if k < n - 1: M[k, k] = -a M[k, k + 1] = a B = Matrix(basis) C = Matrix([[1] + [0]*(n - 1)]) derivs = [eye(n)] for k in xrange(n): derivs.append(M*derivs[k]) l = poly.all_coeffs() l.reverse() res = [0]*n for k, c in enumerate(l): for r, d in enumerate(C*derivs[k]): res[r] += c*d for k, c in enumerate(res): M[n - 1, k] = -c/derivs[n - 1][0, n - 1]/poly.all_coeffs()[0] return Formula(func, z, None, [], B, C, M) else: # Since there are no `ap`, none of the `bq` can be non-positive # integers. basis = [] bq = list(func.bq[:]) for i in range(len(bq)): basis += [hyper([], bq, z)] bq[i] += 1 basis += [hyper([], bq, z)] B = Matrix(basis) n = len(B) C = Matrix([[1] + [0]*(n - 1)]) M = zeros(n) M[0, n - 1] = z/Mul(*func.bq) for k in range(1, n): M[k, k - 1] = func.bq[k - 1] M[k, k] = -func.bq[k - 1] return Formula(func, z, None, [], B, C, M)
def _form_fr(self, fl): """Form the generalized active force. Computes the vector of the generalized active force vector. Used to compute E.o.M. in the form Fr + Fr* = 0. Parameters ========== fl : list Takes in a list of (Point, Vector) or (ReferenceFrame, Vector) tuples which represent the force at a point or torque on a frame. """ if not hasattr(fl, '__iter__'): raise TypeError('Force pairs must be supplied in an iterable.') N = self._inertial self._forcelist = fl[:] u = self._u o = len(u) # number of gen. speeds b = len(fl) # number of forces FR = zeros(o, 1) # pull out relevant velocities for constructing partial velocities vel_list = [] f_list = [] for i in fl: if isinstance(i[0], ReferenceFrame): vel_list += [i[0].ang_vel_in(N)] elif isinstance(i[0], Point): vel_list += [i[0].vel(N)] else: raise TypeError('First entry in pair must be point or frame.') f_list += [i[1]] partials = self._partial_velocity(vel_list, u, N) # Fill Fr with dot product of partial velocities and forces for i in range(o): for j in range(b): FR[i] += partials[j][i] & f_list[j] # In case there are dependent speeds m = len(self._udep) # number of dependent speeds if m != 0: p = o - m FRtilde = FR[:p, 0] FRold = FR[p:o, 0] FRtilde += self._Ars.T * FRold FR = FRtilde self._fr = FR return FR
def balance_reaction(self): """The balance_reaction method will generate a sympy matrix with rows that correspond to each element in the reaction and columns that correspond to each compound. Using the sympy linsolve function it will begin to solve the matrix for each coefficient in terms of of of the compounds. This compound can be assumed to equal 1 in the next method which will allow for scaling of the coefficients while keeping the reaction balanced.""" # The rows of the reaction matrix is each element in the reaction # The columns are the total compounds present in the equation # The equation matrix is defined using the sympy zeros function, all slots are set to zero for easy manipulation reaction_matrix = sp.zeros(len(self.elements), len(self.species)) # The symbols_list contains every compound in the reaction as a sympy symbol so that it can solve the matrix symbols_list = [] # Every compound in the reaction is checked for elements for compound in self.species.keys(): # The periodictable method atoms is called on the compound formula # This identifies how many of each element are present in the compound composition = ptable.formula(compound).atoms # The compound is made into a symbol and stored in the symbols_list symbols_list.append(sp.sympify(compound)) # Every element in the reaction must be considered, even if they are not all present in the compound for element in self.elements: # if the element is not present in the compound, it will be stored as a zero in the composition if element not in composition: composition[element] = 0 # The reaction matrix for the compound and element is set to the compound amount of that element # self.elements.index(element) returns an integer number for the element position in the matrix # list(self.species.keys()).index(compound) returns an integer number for the compound position reaction_matrix[self.elements.index(element), list(self.species.keys()).index(compound)] = composition[element] # Sympy converts the equation matrix to a sympy compatible matrix object # It also creates the solution vector of all 0 reaction_matrix = sp.Matrix(reaction_matrix) solution_vector = sp.Matrix(sp.zeros(len(self.elements), 1)) # A sympy set object containing values is output by the linsolve function from sympy, # which solves the matrix in terms of one of the compounds present (usually the last compound) solution_set = sp.linsolve((reaction_matrix, solution_vector), symbols_list) # This converts the solution_set to usable coefficients self.interpret_balance(symbols_list, solution_set)
def hermitian_basis(dim): """ Returns a basis of the hermitian matrices of size ``dim`` that is orthogonal w.r.t. the Frobenius scalar product. :param dim: size of the matrices :type dim: int Example: >>> import kdotp_symmetry as kp >>> kp.hermitian_basis(2) [Matrix([ [1, 0], [0, 0]]), Matrix([ [0, 0], [0, 1]]), Matrix([ [0, 1], [1, 0]]), Matrix([ [0, -I], [I, 0]])] """ basis = [] # diagonal entries for i in range(dim): mat = sp.zeros(dim) mat[i, i] = 1 basis.append(mat) # off-diagonal entries for i in range(dim): for j in range(i + 1, dim): # real mat = sp.zeros(dim) mat[i, j] = 1 mat[j, i] = 1 basis.append(mat) # imag mat = sp.zeros(dim) mat[i, j] = -sp.numbers.I # pylint: disable=no-member mat[j, i] = sp.numbers.I # pylint: disable=no-member basis.append(mat) # check ONB property assert len(basis) == dim**2 _assert_orthogonal(basis) return basis
def build_hypergeometric_formula(func): """ Create a formula object representing the hypergeometric function ``func``. """ # We know that no `ap` are negative integers, otherwise "detect poly" # would have kicked in. However, `ap` could be empty. In this case we can # use a different basis. # I'm not aware of a basis that works in all cases. from sympy import zeros, Matrix, eye z = Dummy('z') if func.ap: afactors = [_x + a for a in func.ap] bfactors = [_x + b - 1 for b in func.bq] expr = _x*Mul(*bfactors) - z*Mul(*afactors) poly = Poly(expr, _x) n = poly.degree() basis = [] M = zeros(n) for k in range(n): a = func.ap[0] + k basis += [hyper([a] + list(func.ap[1:]), func.bq, z)] if k < n - 1: M[k, k] = -a M[k, k + 1] = a B = Matrix(basis) C = Matrix([[1] + [0]*(n - 1)]) derivs = [eye(n)] for k in range(n): derivs.append(M*derivs[k]) l = poly.all_coeffs() l.reverse() res = [0]*n for k, c in enumerate(l): for r, d in enumerate(C*derivs[k]): res[r] += c*d for k, c in enumerate(res): M[n - 1, k] = -c/derivs[n - 1][0, n - 1]/poly.all_coeffs()[0] return Formula(func, z, None, [], B, C, M) else: # Since there are no `ap`, none of the `bq` can be non-positive # integers. basis = [] bq = list(func.bq[:]) for i in range(len(bq)): basis += [hyper([], bq, z)] bq[i] += 1 basis += [hyper([], bq, z)] B = Matrix(basis) n = len(B) C = Matrix([[1] + [0]*(n - 1)]) M = zeros(n) M[0, n - 1] = z/Mul(*func.bq) for k in range(1, n): M[k, k - 1] = func.bq[k - 1] M[k, k] = -func.bq[k - 1] return Formula(func, z, None, [], B, C, M)
def split_by_ems(self, same_react = False, warn = False): """Split reactions according to the elementary modes they take part in. Return a list of reactions grouped by elementary mode, and the list of reactions not taking part in any elementary mode. :param same_react: if True, do not split reactions with the same reactant. :type same_react: boolean :param warn: if True, give a warning if not all reactions take part in at least one elementary mode. :type warn: boolean :rtype: (list of lists reactions, list of reactions). """ subnet = {} tinvs = self.t_invariants if tinvs.rows == 0: return [self.reactions], [] # case of reactions not taking part in any em if warn: if any(sum(tinvs[:, c]) == 0 for c in range(tinvs.cols)): warnings.warn("Network not covered by elementary modes.") a = sp.SparseMatrix(sp.zeros(self.n_reactions)) for t in range(tinvs.rows): inds = [r for r in range(self.n_reactions) if tinvs[t, r] != 0] for i in inds: for j in inds: a[i, j] = 1 if same_react: for c in self.complexes: inds = [r for r in range(self.n_reactions) if self.reactions[r].reactant == c] for i in inds: for j in inds: a[i, j] = 1 ncc, cc = connected_components(csr_matrix(np.array(a.tolist()).astype(np.int))) rcc = [[self.reactions[r] for r in range(self.n_reactions) if cc[r] == l] for l in range(ncc)] return [rc for rc in rcc if len(rc) > 1], [rc[0] for rc in rcc if len(rc) == 1] ### Print ###
