Python sympy.denom函数代码示例

def tustin(tf, sRate):
	s =	getMapping(tf)['s']
	T = 1/sRate
	poles = sympy.roots(sympy.denom(tf), s, multiple=True)
	centroid = np.mean(np.abs(poles))
	invz = sympy.Symbol('invz', real=True)
	# normalized center frequency derived from poles of filter SISO transfer function
	w0 = 2*np.pi*centroid/(sRate*2*np.pi)
	# modified bilinear transform w/ frequency warping
	bt = w0/sympy.tan(w0*T/2) * ((1-invz)/(1+invz))
	dt = sympy.simplify(tf.subs({s: bt}))
	b = sympy.Poly(sympy.numer(dt)).all_coeffs()[::-1]
	a = sympy.Poly(sympy.denom(dt)).all_coeffs()[::-1]
	normalize = lambda x: float(x/a[0])
	return (map(normalize, b), map(normalize, a))

def poles(G):
    '''
    Return the poles of a multivariable transfer function system. Applies
    Theorem 4.4 (p135).

    Parameters
    ----------
    G : numpy matrix (n x n)
        The transfer function G(s) of the system.

    Returns
    -------
    zero : array
        List of zeros.

    Example
    -------
    >>> def G(s):
    ...     return 1 / (s + 2) * numpy.matrix([[s - 1,  4],
    ...                                       [4.5, 2 * (s - 1)]])
    >>> poles(G)
    [-2.00000000000000]

    Note
    ----
    Not applicable for a non-squared plant, yet.
    '''

    s = sympy.Symbol('s')
    G = sympy.Matrix(G(s))  # convert to sympy matrix object
    det = sympy.simplify(G.det())
    pole = sympy.solve(sympy.denom(det))
    return pole

def find_simple_recurrence_vector(l):
    """
    This function is used internally by other functions from the
    sympy.concrete.guess module. While most users may want to rather use the
    function find_simple_recurrence when looking for recurrence relations
    among rational numbers, the current function may still be useful when
    some post-processing has to be done.

    The function returns a vector of length n when a recurrence relation of
    order n is detected in the sequence of rational numbers v.

    If the returned vector has a length 1, then the returned value is always
    the list [0], which means that no relation has been found.

    While the functions is intended to be used with rational numbers, it should
    work for other kinds of real numbers except for some cases involving
    quadratic numbers; for that reason it should be used with some caution when
    the argument is not a list of rational numbers.

    Examples
    ========

    >>> from sympy.concrete.guess import find_simple_recurrence_vector
    >>> from sympy import fibonacci
    >>> find_simple_recurrence_vector([fibonacci(k) for k in range(12)])
    [1, -1, -1]

    See also
    ========

    See the function sympy.concrete.guess.find_simple_recurrence which is more
    user-friendly.

    """
    q1 = [0]
    q2 = [Integer(1)]
    b, z = 0, len(l) >> 1
    while len(q2) <= z:
        while l[b]==0:
            b += 1
            if b == len(l):
                c = 1
                for x in q2:
                    c = lcm(c, denom(x))
                if q2[0]*c < 0: c = -c
                for k in range(len(q2)):
                    q2[k] = int(q2[k]*c)
                return q2
        a = Integer(1)/l[b]
        m = [a]
        for k in range(b+1, len(l)):
            m.append(-sum(l[j+1]*m[b-j-1] for j in range(b, k))*a)
        l, m = m, [0] * max(len(q2), b+len(q1))
        for k in range(len(q2)):
            m[k] = a*q2[k]
        for k in range(b, b+len(q1)):
            m[k] += q1[k-b]
        while m[-1]==0: m.pop() # because trailing zeros can occur
        q1, q2, b = q2, m, 1
    return [0]

def e2nd(expression):
	""" basic helper function that accepts a sympy expression, expands it,
	attempts to simplify it, and returns a numerator and denomenator pair for the instantiation of a scipy
	LTI system object. """
	expression = expression.expand()
	expression = expression.cancel()
	n = sympy.Poly(sympy.numer(expression), s).all_coeffs()
	d = sympy.Poly(sympy.denom(expression), s).all_coeffs()
	n = map(float, n)
	d = map(float, d)
	return (n, d)

    def __new__(cls, b, e, evaluate=True):
        from sympy.functions.elementary.exponential import exp_polar
        from sympy.functions import log
        # don't optimize "if e==0; return 1" here; it's better to handle that
        # in the calling routine so this doesn't get called
        b = _sympify(b)
        e = _sympify(e)
        if evaluate:
            if e is S.Zero:
                return S.One
            elif e is S.One:
                return b
            elif S.NaN in (b, e):
                if b is S.One:  # already handled e == 0 above
                    return S.One
                return S.NaN
            else:
                if e.func == log:
                    if len(e.args) == 2:
                        lbase = e.args[1]
                    else:
                        lbase = S.Exp1
                    if lbase == b:
                        return e.args[0]

                if e is Mul and e.args[1].func == log:
                    if len(e.args[1].args) == 2:
                        lbase = e.args[1].args[1]
                    else:
                        lbase = S.Exp1
                    if lbase == b:
                        return e.args[1].args[0] ** e.args[0]

                # recognize base as E
                if not e.is_Atom and b is not S.Exp1 and b.func is not exp_polar:
                    from sympy import numer, denom, log, sign, im, factor_terms
                    c, ex = factor_terms(e, sign=False).as_coeff_Mul()
                    den = denom(ex)
                    if den.func is log and den.args[0] == b:
                        return S.Exp1**(c*numer(ex))
                    elif den.is_Add:
                        s = sign(im(b))
                        if s.is_Number and s and den == \
                                log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
                            return S.Exp1**(c*numer(ex))

                obj = b._eval_power(e)
                if obj is not None:
                    return obj
        obj = Expr.__new__(cls, b, e)
        obj.is_commutative = (b.is_commutative and e.is_commutative)
        return obj

    def Pow(expr, assumptions):
        """
        Imaginary**integer/odd  -> Imaginary
        Imaginary**integer/even -> Real if integer % 2 == 0
        b**Imaginary            -> !Imaginary if exponent is an integer multiple of I*pi/log(b)
        Imaginary**Real         -> ?
        Negative**even root     -> Imaginary
        Negative**odd root      -> Real
        Negative**Real          -> Imaginary
        Real**Integer           -> Real
        Real**Positive          -> Real
        """
        if expr.is_number:
            return AskImaginaryHandler._number(expr, assumptions)

        if expr.base.func == C.exp:
            if ask(Q.imaginary(expr.base.args[0]), assumptions):
                if ask(Q.imaginary(expr.exp), assumptions):
                    return False
                i = expr.base.args[0]/I/pi
                if ask(Q.integer(2*i), assumptions):
                    return ask(Q.imaginary(((-1)**i)**expr.exp), assumptions)

        if ask(Q.imaginary(expr.base), assumptions):
            if ask(Q.integer(expr.exp), assumptions):
                odd = ask(Q.odd(expr.exp), assumptions)
                if odd is not None:
                    return odd
                return

        if ask(Q.imaginary(expr.exp), assumptions):
            imlog = ask(Q.imaginary(C.log(expr.base)), assumptions)
            if imlog is not None:
                return False  # I**i -> real; (2*I)**i -> complex ==> not imaginary

        if ask(Q.real(expr.base), assumptions):
            if ask(Q.real(expr.exp), assumptions):
                if ask(Q.rational(expr.exp) & Q.even(denom(expr.exp)), assumptions):
                    return ask(Q.negative(expr.base), assumptions)
                elif ask(Q.integer(expr.exp), assumptions):
                    return False
                elif ask(Q.positive(expr.base), assumptions):
                    return False
                elif ask(Q.negative(expr.base), assumptions):
                    return True

def mutate_rat(value, *, keep_sign=True, mutate_up=None):
    '''Mutate rational number'''

    r = random.random()
    numer, denom = sp.numer(value), sp.denom(value)

    # swap numerator and denominator
    if r < 0.20:
        numer, denom = (mutate_int(denom, mutate_up=mutate_up, keep_sign=keep_sign),
                        mutate_int(numer, mutate_up=mutate_up))
    # mutate numerator
    elif r < 0.65:
        numer = mutate_int(numer, mutate_up=mutate_up, keep_sign=keep_sign)
    # mutate denominator
    else:
        denom = mutate_int(denom, mutate_up=mutate_up, keep_sign=keep_sign)

    return sp.Rational(numer, denom or 1)

    def __new__(cls, b, e, evaluate=None):
        if evaluate is None:
            evaluate = global_evaluate[0]
        from sympy.functions.elementary.exponential import exp_polar

        b = _sympify(b)
        e = _sympify(e)
        if evaluate:
            if e is S.Zero:
                return S.One
            elif e is S.One:
                return b
            elif e.is_integer and _coeff_isneg(b):
                if e.is_even:
                    b = -b
                elif e.is_odd:
                    return -Pow(-b, e)
            if b is S.One:
                if e in (S.NaN, S.Infinity, -S.Infinity):
                    return S.NaN
                return S.One
            elif S.NaN in (b, e):  # XXX S.NaN**x -> S.NaN under assumption that x != 0
                return S.NaN
            else:
                # recognize base as E
                if not e.is_Atom and b is not S.Exp1 and b.func is not exp_polar:
                    from sympy import numer, denom, log, sign, im, factor_terms
                    c, ex = factor_terms(e, sign=False).as_coeff_Mul()
                    den = denom(ex)
                    if den.func is log and den.args[0] == b:
                        return S.Exp1**(c*numer(ex))
                    elif den.is_Add:
                        s = sign(im(b))
                        if s.is_Number and s and den == \
                                log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
                            return S.Exp1**(c*numer(ex))

                obj = b._eval_power(e)
                if obj is not None:
                    return obj
        obj = Expr.__new__(cls, b, e)
        obj.is_commutative = (b.is_commutative and e.is_commutative)
        return obj

def symmetryDetection(allVariables, diffEquations, observables, obsFunctions, initFunctions,
						predictions, predFunctions, ansatz = 'uni', pMax = 2, inputs = [], 
						fixed = [], parallel = 1, allTrafos = False, timeTrans = False,
						pretty = True, suffix=''):
	
	n = len(allVariables)
	m = len(diffEquations)
	h = len(observables)

	###########################################################################################
	#############################     prepare equations    ####################################
	###########################################################################################
	sys.stdout.write('Preparing equations...')
	sys.stdout.flush()

	# make infinitesimal ansatz
	infis, diffInfis, rs = makeAnsatz(ansatz, allVariables, m, len(inputs), pMax, fixed)

	# get infinitesimals of time transformation
	if timeTrans:
		rs.append(spy.var('r_T_1'))
		diffInfiT = rs[-1]
		allVariables += [T]
	else:
		diffInfiT = None

	# and convert to polynomial
	infis, diffInfis = transformInfisToPoly(infis, diffInfis, allVariables, rs, parallel, ansatz)
	
	diffInfiT = Apoly(diffInfiT, allVariables, rs)

	### extract numerator and denominator of equations
	#differential equations
	numerators = [0]*m
	denominators = [0]*m
	for k in range(m):
		rational = spy.together(diffEquations[k])
		numerators[k] = Apoly(spy.numer(rational), allVariables, None)
		denominators[k] = Apoly(spy.denom(rational), allVariables, None)

	#observation functions
	obsNumerators = [0]*h
	obsDenominatros = [0]*h
	for k in range(h):
		rational = spy.together(obsFunctions[k])
		obsNumerators[k] = Apoly(spy.numer(rational), allVariables, None)
		obsDenominatros[k] = Apoly(spy.denom(rational), allVariables, None)

	#initial functions
	if len(initFunctions) != 0:
		initNumerators = [0]*m
		initDenominatros = [0]*m
		for k in range(m):
			rational = spy.together(initFunctions[k])
			initNumerators[k] = Apoly(spy.numer(rational), allVariables, None)
			initDenominatros[k] = Apoly(spy.denom(rational), allVariables, None)
	else:
		initNumerators = []
		initDenominatros = []

	### calculate numerator of derivatives of equations
	#differential equatioins
	derivativesNum = [0]*m
	for i in range(m):
		derivativesNum[i] = [0]*n
	for k in range(m):
		for l in range(n):
			derivativesNum[k][l] = Apoly(None, allVariables, None)
			derivativesNum[k][l].add(numerators[k].diff(l).mul(denominators[k]))
			derivativesNum[k][l].sub(numerators[k].mul(denominators[k].diff(l)))

	#observation functions
	obsDerivativesNum = [0]*h
	for i in range(h):
		obsDerivativesNum[i] = [0]*n
	for k in range(h):
		for l in range(n):
			obsDerivativesNum[k][l] = Apoly(None, allVariables, None)
			obsDerivativesNum[k][l].add(obsNumerators[k].diff(l).mul(obsDenominatros[k]))
			obsDerivativesNum[k][l].sub(obsNumerators[k].mul(obsDenominatros[k].diff(l)))

	#initial functions
	if len(initFunctions) != 0:
		initDerivativesNum = [0]*len(initFunctions)
		for i in range(m):
			initDerivativesNum[i] = [0]*n
		for k in range(m):
			for l in range(n):
				initDerivativesNum[k][l] = Apoly(None, allVariables, None)
				initDerivativesNum[k][l].add(initNumerators[k].diff(l).mul(initDenominatros[k]))
				initDerivativesNum[k][l].sub(initNumerators[k].mul(initDenominatros[k].diff(l)))
	else:
		initDerivativesNum = []

	sys.stdout.write('\rPreparing equations...done\n')
	sys.stdout.flush()

	###########################################################################################
	############################     build linear system    ###################################
	###########################################################################################
#.........这里部分代码省略.........

m2_2 = G[:,[0,2]]
m2 = m2_1.row_join(m2_2)
m2_3 = G[:,[0,1]]
m2 = m2.row_join(m2_3)
print(m2_1)
print(m2_2)
print(m2_3)

# combine order one and two minors

print('---All poles---')
m1count = np.shape(m1)[0]
#find out how many roots there are in total between all 1st order minors
n1Roots = 0
for m in range(m1count):
    n1Roots = n1Roots + sp.degree(sp.denom(m1[m]),s)

#find out how many roots there are in total between all 2nd order minors
n2Roots = 0
r,c = np.shape(m2)
for i in range(r):
    for j in range(c):
        n2Roots = n2Roots + sp.degree(sp.denom(m2[i,j]),s)
print(n2Roots)

# calculate and find common roots 
roots_denom = [None]*(n1Roots + n2Roots)
n = 0
for m in range(m1count):
    denom_roots = sp.solve(sp.denom(m1[m]))
    for i  in range(len(denom_roots)):

def test_issue_5933():
    from sympy import Polygon, RegularPolygon, denom
    x = Polygon(*RegularPolygon((0, 0), 1, 5).vertices).centroid.x
    assert abs(denom(x).n()) > 1e-12
    assert abs(denom(radsimp(x))) > 1e-12  # in case simplify didn't handle it

def gen_lobatto(max_order):
    assert max_order > 2

    x = sm.symbols('x')

    lobs = [0, 1]
    lobs[0] = (1 - x) / 2
    lobs[1] = (1 + x) / 2

    dlobs = [lob.diff('x') for lob in lobs]

    legs = [sm.legendre(0, 'y')]
    clegs = [sm.ccode(legs[0])]
    dlegs = [sm.legendre(0, 'y').diff('y')]
    cdlegs = [sm.ccode(dlegs[0])]

    clobs = [sm.ccode(lob) for lob in lobs]
    cdlobs = [sm.ccode(dlob) for dlob in dlobs]

    denoms = [] # for lobs.

    for ii in range(2, max_order + 1):
        coef = sm.sympify('sqrt(2 * (2 * %s - 1)) / 2' % ii)
        leg = sm.legendre(ii - 1, 'y')

        pleg = leg.as_poly()
        coefs = pleg.all_coeffs()
        denom = max(sm.denom(val) for val in coefs)

        cleg = sm.ccode(sm.horner(leg*denom)/denom)

        dleg = leg.diff('y')
        cdleg = sm.ccode(sm.horner(dleg*denom)/denom)

        lob = sm.simplify(coef * sm.integrate(leg, ('y', -1, x)))
        lobnc = sm.simplify(sm.integrate(leg, ('y', -1, x)))

        plobnc = lobnc.as_poly()
        coefs = plobnc.all_coeffs()
        denom = sm.denom(coef) * max(sm.denom(val) for val in coefs)

        clob = sm.ccode(sm.horner(lob*denom)/denom)

        dlob = lob.diff('x')
        cdlob = sm.ccode(sm.horner(dlob*denom)/denom)

        legs.append(leg)
        clegs.append(cleg)
        dlegs.append(dleg)
        cdlegs.append(cdleg)
        lobs.append(lob)
        clobs.append(clob)
        dlobs.append(dlob)
        cdlobs.append(cdlob)
        denoms.append(denom)

    coef = sm.sympify('sqrt(2 * (2 * %s - 1)) / 2' % (max_order + 1))
    leg = sm.legendre(max_order, 'y')

    pleg = leg.as_poly()
    coefs = pleg.all_coeffs()
    denom = max(sm.denom(val) for val in coefs)

    cleg = sm.ccode(sm.horner(leg*denom)/denom)

    dleg = leg.diff('y')
    cdleg = sm.ccode(sm.horner(dleg*denom)/denom)

    legs.append(leg)
    clegs.append(cleg)
    dlegs.append(dleg)
    cdlegs.append(cdleg)

    kerns = []
    ckerns = []
    dkerns = []
    cdkerns = []
    for ii, lob in enumerate(lobs[2:]):
        kern = sm.simplify(lob / (lobs[0] * lobs[1]))
        dkern = kern.diff('x')

        denom = denoms[ii] / 4
        ckern = sm.ccode(sm.horner(kern*denom)/denom)
        cdkern = sm.ccode(sm.horner(dkern*denom)/denom)

        kerns.append(kern)
        ckerns.append(ckern)
        dkerns.append(dkern)
        cdkerns.append(cdkern)

    return (legs, clegs, dlegs, cdlegs,
            lobs, clobs, dlobs, cdlobs,
            kerns, ckerns, dkerns, cdkerns,
            denoms)

	def simplify(self):
		self._tuple = tuple(sympy.factor(ex) for ex in self._tuple)
		m = functools.reduce(lambda x, y: x * sympy.denom(y), self._tuple, 1)
		self._tuple = tuple(ex * m for ex in self._tuple)
		d = functools.reduce(sympy.gcd, self._tuple)
		self._tuple = tuple(sympy.factor(ex / d) for ex in self._tuple)

tf = ahkabHelpers.reduceTF(tf, mycircuit)

sRate = 44.1e3

fs = lambda x: sympy.N(abs(tf.subs({s:sympy.I*x})))
ws = np.logspace(1, np.log10(sRate), 5e2)
mags = map(fs,ws)

print tf
(b, a) = ahkabHelpers.tustin(tf, sRate)
print b,a

pylab.semilogx(ws, map(fs, ws), 'v', label="from transfer function")
pylab.semilogx(r['ac']['w'][::10], np.abs(r['ac']['VU1o'][::10]), '-', label='from AC simulation')
pylab.vlines(np.abs(sympy.roots(sympy.denom(tf), s, multiple=True)), 0, 1, 'r')

# build white noise input vector, normalized to \pm 1.0
x = list(2*np.random.random(sRate)-1.0)
# allocate output vector for y[n-2] indexing to work
y = [0.0]*len(x)

# Direct Form I
for n in range(len(x)):
	y[n] = b[0]*x[n] + b[1]*x[n-1] + b[2]*x[n-2] - a[1]*y[n-1] - a[2]*y[n-2]

# frequencies to radians
freqs = (np.fft.fftfreq(len(y), 1/sRate))[0:len(y)/2] * (2*np.pi)
mags = 	(np.abs(np.fft.fft(np.array(y))))[0:len(y)/2]

# normalize magnitude to 1.0

def ac_analysis(param_d, param_l, instance, file_sufix):
    """ Performs ac analysis

        param_d: substitutions for symbols, or named parameters for plot
        
        param_l: expresions to plot

        format is:
        .ac expresion0 [expresion1 expresion2 ...] sweep = parameter_to_sweep [symmbol_or_option0 = value0
        symmbol_or_option1 = value1 ...]
        
        expresions are on positiona parameters list can hold expresions containing parameters or functions of nodal
        voltages [v(node)] and element and port currents [i(element) isub(port)]

        named parameters (param_d) can contain options for analysis or substitutions for symbols, substituions will be
        done as they are without any parsing, be aware of symbol and option names clashes.

        Config options:
        fstart:         first value of frequency for AC analysis [float]
        fstop:          last value of frequency for AC analysis [float]
        fscale:         scale for frequency points [linear | log]
        npoints:        numbers of points for ac analysis [integer]
        yscale:         scale for y-axis to being displayed on [linear or log]
        hold:           hold plot for next analysis and don't save it to file [yes | no]
        show_poles:     show poles of function on plot [yes | no]
        show_zeroes:    show zeros of function on plot [yes | no]
        title:          display title above ac plot [string]
        show_legend:    show legend on plot [yes | no]
        xkcd:           style plot to be xkcd like scetch
    """
    warnings.filterwarnings('ignore')  # Just getting rid of those fake casting from complex warnings
    s, w = sympy.symbols(('s', 'w'))
    config = {'fstart': 1,
              'fstop': 1e6,
              'fscale': 'log',
              'npoints': 100,
              'yscale': 'log',
              'type': 'amp',
              'hold': 'no',
              'show_poles': 'yes',
              'show_zeros': 'yes',
              'title': None,
              'show_legend': 'no',
              'xkcd': 'no'}

    for config_name in config.keys():
        if config_name in param_d:
            config.update({config_name: param_d[config_name]})
            param_d.pop(config_name)

    subst = []
    for symbol, value in param_d.iteritems():
        tokens = scs_parser.parse_param_expresion(value)
        try:
            value = float(sympy.sympify(scs_parser.params2values(tokens, instance.paramsd),sympy.abc._clash))
        except ValueError:
            raise scs_errors.ScsAnalysisError("Passed subsitution for %s is not a number")

        subst.append((symbol, value))

    if config['fscale'] == 'log':
        fs = np.logspace(np.log10(float(config['fstart'])), np.log10(float(config['fstop'])), int(config['npoints']))
    elif config['fscale'] == 'linear':
        fs = np.linspace(float(config['fstart']), float(config['fstop']), int(config['npoints']))
    else:
        raise scs_errors.ScsAnalysisError(("Option %s for fscale invalid!" % config['yscale']))

    if config['yscale'] != 'log' and config['yscale'] != 'linear':
        raise scs_errors.ScsAnalysisError(("Option %s for yscale invalid!" % config['fscale']))

    filename = "%s.results" % file_sufix

    with open(filename, 'a') as fil:
        if config['xkcd'] == 'yes':
            plt.xkcd()
        plt.hold(True)

        for expresion in param_l:
            fil.write("%s: %s \n---------------------\n" % ('AC analysis of', expresion))
            tokens = scs_parser.parse_analysis_expresion(expresion)
            value0 = sympy.factor(sympy.sympify(scs_parser.results2values(tokens, instance),sympy.abc._clash), s).simplify()
            fil.write("%s = %s \n\n" % (expresion, str(value0)))
            denominator = sympy.denom(value0)
            numerator = sympy.numer(value0)
            poles = sympy.solve(denominator, s)
            zeros = sympy.solve(numerator, s)
            poles_r = sympy.roots(denominator, s)
            zeros_r = sympy.roots(numerator, s)
            gdc = str(value0.subs(s, 0).simplify())
            fil.write('G_DC = %s\n\n' % gdc)

            p = 0
            titled = 1
            for pole, degree in poles_r.iteritems():
                if pole == 0:
                    titled *= s ** degree
                else:
                    titled *= (s / sympy.symbols("\\omega_p%d" % p) + 1)
                p += 1
            z = 0
#.........这里部分代码省略.........