From da0919c520f6f4595718599efdca1722d63dbf5d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner <inform@tiker.net> Date: Mon, 13 Nov 2006 20:20:04 -0500 Subject: [PATCH] [pymbolic @ inform@tiker.net-20061114012004-d027ae58069d69bd] First test already running again after Constant-ripout, but differentiation still buggy. --- setup.py | 2 +- src/__init__.py | 5 +- src/compiler.py | 5 +- src/mapper/__init__.py | 94 +++++++++++++++++++++--------------- src/mapper/dependency.py | 7 ++- src/mapper/differentiator.py | 18 +++---- src/mapper/evaluator.py | 2 +- src/mapper/stringifier.py | 24 +++++++-- src/parser.py | 6 +-- src/polynomial.py | 35 +++----------- src/primitives.py | 66 +++++++++++-------------- src/rational.py | 15 +++++- src/traits.py | 5 +- 13 files changed, 147 insertions(+), 137 deletions(-) diff --git a/setup.py b/setup.py index beca570..6119830 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ import os.path setup(name="pymbolic", version="0.10", description="A package for symbolic computation", - author=u"Andreas Kl�ckner", + author=u"Andreas Kloeckner", author_email="inform@tiker.net", license = "BSD, like Python itself", url="http://news.tiker.net/software/pymbolic", diff --git a/src/__init__.py b/src/__init__.py index 5ded8ac..6e08062 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -9,7 +9,6 @@ import pymbolic.mapper.differentiator import pymbolic.primitives var = pymbolic.primitives.Variable -const = pymbolic.primitives.Constant sum = pymbolic.primitives.sum subscript = pymbolic.primitives.subscript product = pymbolic.primitives.product @@ -38,7 +37,7 @@ def jacobian(expression_list, variables): return [grad(expr, variables) for expr in expression_list] def laplace(expression, variables): - return sum(*[differentiate(differentiate(expression,var), var) for var in variables]) + return sum(differentiate(differentiate(expression,var), var) for var in variables) @@ -75,7 +74,7 @@ if __name__ == "__main__": print ex print repr(parse("x+y")) print evaluate(ex, {"alpha":5, "math":math, "x":-math.pi}) - compiled = compile(substitute(ex, {var("alpha"): const(5)})) + compiled = compile(substitute(ex, {var("alpha"): 5})) print compiled(-math.pi) import cPickle as pickle pickle.dumps(compiled) diff --git a/src/compiler.py b/src/compiler.py index 349d517..74c4b4b 100644 --- a/src/compiler.py +++ b/src/compiler.py @@ -1,4 +1,3 @@ -import sets import math import pymbolic @@ -22,8 +21,8 @@ class CompiledExpression: ctx = self.context() used_variables = pymbolic.get_dependencies(self._Expression) - used_variables -= sets.Set(self._Variables) - used_variables -= sets.Set(pymbolic.var(key) for key in ctx.keys()) + used_variables -= set(self._Variables) + used_variables -= set(pymbolic.var(key) for key in ctx.keys()) used_variables = list(used_variables) used_variables.sort() all_variables = self._Variables + used_variables diff --git a/src/mapper/__init__.py b/src/mapper/__init__.py index 2f71ad6..a2ee13e 100644 --- a/src/mapper/__init__.py +++ b/src/mapper/__init__.py @@ -1,12 +1,12 @@ class Mapper(object): def __call__(self, victim, *args, **kwargs): try: - victim.invoke_mapper(self, *args, **kwargs) + return victim.invoke_mapper(self, *args, **kwargs) except AttributeError: - self.map_constant(victim) + return self.map_constant(victim, *args, **kwargs) - def map_rational(self, expr): - return self.map_quotient(self, expr) + def map_rational(self, expr, *args, **kwargs): + return self.map_quotient(expr, *args, **kwargs) @@ -16,30 +16,39 @@ class CombineMapper(Mapper): def combine(self, values): raise NotImplementedError - def map_call(self, expr): - return self.combine([self(expr.function)] + - [self(child) for child in expr.parameters]) + def map_call(self, expr, *args, **kwargs): + return self.combine( + (self(expr.function, *args, **kwargs),) + + tuple( + self(child, *args, **kwargs) for child in expr.parameters) + ) - def map_subscript(self, expr): + def map_subscript(self, expr, *args, **kwargs): return self.combine( - [self(expr.aggregate), self(expr.index)]) + [self(expr.aggregate, *args, **kwargs), + self(expr.index, *args, **kwargs)]) - def map_lookup(self, expr): - return self(expr.aggregate) + def map_lookup(self, expr, *args, **kwargs): + return self(expr.aggregate, *args, **kwargs) - def map_negation(self, expr): - return self(expr.child) + def map_negation(self, expr, *args, **kwargs): + return self(expr.child, *args, **kwargs) - def map_sum(self, expr): - return self.combine(self(child) for child in expr.children) + def map_sum(self, expr, *args, **kwargs): + return self.combine(self(child, *args, **kwargs) + for child in expr.children) map_product = map_sum - def map_quotient(self, expr): - return self.combine((self(expr.numerator), self(expr.denominator))) + def map_quotient(self, expr, *args, **kwargs): + return self.combine(( + self(expr.numerator, *args, **kwargs), + self(expr.denominator, *args, **kwargs))) - def map_power(self, expr): - return self.combine(self(expr.base), self(expr.exponent)) + def map_power(self, expr, *args, **kwargs): + return self.combine(( + self(expr.base, *args, **kwargs), + self(expr.exponent, *args, **kwargs))) map_list = map_sum @@ -48,42 +57,49 @@ class CombineMapper(Mapper): class IdentityMapper(Mapper): - def map_constant(self, expr): + def map_constant(self, expr, *args, **kwargs): + # leaf -- no need to rebuild return expr - def map_variable(self, expr): + def map_variable(self, expr, *args, **kwargs): + # leaf -- no need to rebuild return expr - def map_call(self, expr): + def map_call(self, expr, *args, **kwargs): return expr.__class__( - self(expr.function), - tuple(self(child) + self(expr.function, *args, **kwargs), + tuple(self(child, *args, **kwargs) for child in expr.parameters)) - def map_subscript(self, expr): - return expr.__class__(self(expr.aggregate), self(expr.index)) + def map_subscript(self, expr, *args, **kwargs): + return expr.__class__( + self(expr.aggregate, *args, **kwargs), + self(expr.index, *args, **kwargs)) - def map_lookup(self, expr): - return expr.__class__(self(expr.aggregate), expr.name) + def map_lookup(self, expr, *args, **kwargs): + return expr.__class__( + self(expr.aggregate, *args, **kwargs), + expr.name) - def map_negation(self, expr): - return expr.__class__(self(expr.child)) + def map_negation(self, expr, *args, **kwargs): + return expr.__class__(self(expr.child, *args, **kwargs)) - def map_sum(self, expr): + def map_sum(self, expr, *args, **kwargs): return expr.__class__( - *[self(child) for child in expr.children]) + *[self(child, *args, **kwargs) + for child in expr.children]) map_product = map_sum - def map_rational(self, expr): - return expr.__class__(self(expr.numerator), - self(expr.denominator)) + def map_rational(self, expr, *args, **kwargs): + return expr.__class__(self(expr.numerator, *args, **kwargs), + self(expr.denominator, *args, **kwargs)) - def map_power(self, expr): - return expr.__class__(self(expr.base), - self(expr.exponent)) + def map_power(self, expr, *args, **kwargs): + return expr.__class__(self(expr.base, *args, **kwargs), + self(expr.exponent, *args, **kwargs)) - def map_polynomial(self, expr): + def map_polynomial(self, expr, *args, **kwargs): raise NotImplementedError map_list = map_sum diff --git a/src/mapper/dependency.py b/src/mapper/dependency.py index cb84a71..11cdab9 100644 --- a/src/mapper/dependency.py +++ b/src/mapper/dependency.py @@ -1,4 +1,3 @@ -import sets import pymbolic.mapper @@ -10,10 +9,10 @@ class DependencyMapper(pymbolic.mapper.CombineMapper): return reduce(operator.or_, values) def map_constant(self, expr): - return sets.Set() + return set() def map_variable(self, expr): - return sets.Set([expr]) + return set([expr]) @@ -28,5 +27,5 @@ def is_constant(expr, with_respect_to=None): if not with_respect_to: return not bool(get_dependencies(expr)) else: - return not (sets.Set(with_respect_to) and get_dependencies(expr)) + return not (set(with_respect_to) and get_dependencies(expr)) diff --git a/src/mapper/differentiator.py b/src/mapper/differentiator.py index 855cf26..1111a21 100644 --- a/src/mapper/differentiator.py +++ b/src/mapper/differentiator.py @@ -24,7 +24,7 @@ def map_math_functions_by_name(i, func, pars): elif f is math.tan and len(pars) == 1: return make_f("tan")(*pars)**2+1 elif f is math.log and len(pars) == 1: - return primitives.Constant(1)/pars[0] + return primitives.quotient(1, pars[0]) elif f is math.exp and len(pars) == 1: return make_f("exp")(*pars) else: @@ -39,20 +39,20 @@ class DifferentiationMapper(pymbolic.mapper.Mapper): self.FunctionMap = func_map def map_constant(self, expr): - return primitives.Constant(0) + return 0 def map_variable(self, expr): if expr == self.Variable: - return primitives.Constant(1) + return 1 else: - return primitives.Constant(0) + return 0 def map_call(self, expr): - return pymbolic.sum(*( + return pymbolic.sum( self.FunctionMap(i, expr.function, expr.parameters) * self(par) for i, par in enumerate(expr.parameters) - if not self._isc(par))) + if not self._isc(par)) map_subscript = map_variable @@ -79,7 +79,7 @@ class DifferentiationMapper(pymbolic.mapper.Mapper): g_const = self._isc(g) if f_const and g_const: - return primitives.Constant(0) + return 0 elif f_const: return -f*self(g)/g**2 elif g_const: @@ -93,10 +93,10 @@ class DifferentiationMapper(pymbolic.mapper.Mapper): f_const = self._isc(f) g_const = self._isc(g) - log = primitives.Constant("log") + log = "log" if f_const and g_const: - return primitives.Constant(0) + return 0 elif f_const: return log(f) * f**g * self(g) elif g_const: diff --git a/src/mapper/evaluator.py b/src/mapper/evaluator.py index 70c60dd..9ed9e8e 100644 --- a/src/mapper/evaluator.py +++ b/src/mapper/evaluator.py @@ -14,7 +14,7 @@ class EvaluationMapper(pymbolic.mapper.Mapper): self.Context = context def map_constant(self, expr): - return expr.value + return expr def map_variable(self, expr): try: diff --git a/src/mapper/stringifier.py b/src/mapper/stringifier.py index 80e5051..4ed306e 100644 --- a/src/mapper/stringifier.py +++ b/src/mapper/stringifier.py @@ -14,7 +14,7 @@ PREC_NONE = 0 class StringifyMapper(pymbolic.mapper.Mapper): def map_constant(self, expr, enclosing_prec): - return str(expr.value) + return str(expr) def map_variable(self, expr, enclosing_prec): return expr.name @@ -75,7 +75,6 @@ class StringifyMapper(pymbolic.mapper.Mapper): return "(%s)" % result else: return result - map_rational = map_quotient def map_power(self, expr, enclosing_prec): result = "%s**%s" % ( @@ -88,7 +87,26 @@ class StringifyMapper(pymbolic.mapper.Mapper): return result def map_polynomial(self, expr, enclosing_prec): - raise NotImplementedError + sbase = str(expr.base) + def stringify_expcoeff((exp, coeff)): + if exp == 1: + if not (coeff-1): + return sbase + else: + return "%s*%s" % (str(coeff), sbase) + elif exp == 0: + return str(coeff) + else: + if not (coeff-1): + return "%s**%s" % (sbase, exp) + else: + return "%s*%s**%s" % (str(coeff), sbase, exp) + + result = "%s" % " + ".join(stringify_expcoeff(i) for i in expr.data[::-1]) + if enclosing_prec > PREC_SUM: + return "(%s)" % result + else: + return result def map_list(self, expr, enclosing_prec): return "[%s]" % ", ".join([self(i, PREC_NONE) for i in expr.children]) diff --git a/src/parser.py b/src/parser.py index 4231d1b..e733de5 100644 --- a/src/parser.py +++ b/src/parser.py @@ -51,11 +51,11 @@ def parse(expr_str): def parse_terminal(pstate): next_tag = pstate.next_tag() if next_tag is _int: - return primitives.Constant(int(pstate.next_str_and_advance())) + return int(pstate.next_str_and_advance()) elif next_tag is _float: - return primitives.Constant(float(pstate.next_str_and_advance())) + return float(pstate.next_str_and_advance()) elif next_tag is _imaginary: - return primitives.Constant(complex(pstate.next_str_and_advance())) + return complex(pstate.next_str_and_advance()) elif next_tag is _identifier: return primitives.Variable(pstate.next_str_and_advance()) else: diff --git a/src/polynomial.py b/src/polynomial.py index 53495c1..9394418 100644 --- a/src/polynomial.py +++ b/src/polynomial.py @@ -117,8 +117,7 @@ class Polynomial(primitives.Expression): return Polynomial(self.Base, [(exp, coeff * other) for exp, coeff in self.Data]) - result = [] - + result = [] for s_exp, s_coeff in self.Data: for o_exp, o_coeff in other.Data: result.append((s_exp+o_exp, s_coeff*o_coeff)) @@ -169,24 +168,6 @@ class Polynomial(primitives.Expression): def __mod__(self): return self.__divmod__(self, other)[1] - def __str__(self): - sbase = str(self.Base) - def stringify_expcoeff((exp, coeff)): - if exp == 1: - if not (coeff-1): - return sbase - else: - return "%s*%s" % (str(coeff), sbase) - elif exp == 0: - return str(coeff) - else: - if not (coeff-1): - return "%s**%s" % (sbase, exp) - else: - return "%s*%s**%s" % (str(coeff), sbase, exp) - - return "(%s)" % " + ".join(stringify_expcoeff(i) for i in self.Data[::-1]) - def _data(self): return self.Data data = property(_data) @@ -206,10 +187,8 @@ class Polynomial(primitives.Expression): return -1 degree = property(_degree) - def __repr__(self): - return "%s(%s, %s)" % (self.__class__.__name__, - repr(self.Base), - repr(self.Data)) + def __getinitargs__(self): + return (self.Base, self.Data, self.Unit) def invoke_mapper(self, mapper, *args, **kwargs): return mapper.map_polynomial(self, *args, **kwargs) @@ -255,10 +234,9 @@ class PolynomialTraits(traits.EuclideanRingTraits): if __name__ == "__main__": - - x = Polynomial(pymbolic.var("x"), unit=pymbolic.const(1)) - y = Polynomial(pymbolic.var("y"), unit=pymbolic.const(1)) - xpoly = x**2 + pymbolic.const(1) + x = Polynomial(pymbolic.var("x")) + y = Polynomial(pymbolic.var("y")) + xpoly = x**2 + 1 ypoly = -y**2*xpoly + xpoly #print xpoly #print ypoly @@ -270,6 +248,7 @@ if __name__ == "__main__": xp3 = xpoly**3 print xp3 print integrate(xp3) + print differentiate(integrate(xp3)) #print 3*xpoly**3 + 1 #print xpoly diff --git a/src/primitives.py b/src/primitives.py index 56fdc00..551249c 100644 --- a/src/primitives.py +++ b/src/primitives.py @@ -64,13 +64,11 @@ class Expression(object): def __rdiv__(self, other): assert not isinstance(other, Expression) - return quotient(Constant(other), self) + return quotient(other, self) def __pow__(self, other): - if not isinstance(other, Expression): - other = Constant(other) if not other: # exponent zero - return Constant(1) + return 1 elif not (other-1): # exponent one return self return Power(self, other) @@ -78,9 +76,9 @@ class Expression(object): def __rpow__(self, other): assert not isinstance(other, Expression) if not other: # base zero - return Constant(0) + return 0 elif not (other-1): # base one - return Constant(1) + return 1 return Power(other, self) def __neg__(self): @@ -92,7 +90,7 @@ class Expression(object): if isinstance(par, Expression): processed.append(par) else: - processed.append(Constant(par)) + processed.append(par) return Call(self, tuple(processed)) @@ -237,26 +235,20 @@ class Sum(Expression): return isinstance(other, Sum) and (self._Children == other._Children) def __add__(self, other): - if not isinstance(other, Expression): - other = Constant(other) - elif isinstance(other, Sum): + if isinstance(other, Sum): return Sum(*(self._Children + other._Children)) if not other: return self return Sum(*(self._Children + (other,))) def __radd__(self, other): - if not isinstance(other, Expression): - other = Constant(other) - elif isinstance(other, Sum): + if isinstance(other, Sum): return Sum(*(other._Children + self._Children)) if not other: return self return Sum(*((other,) + self._Children)) def __sub__(self, other): - if not isinstance(other, Expression): - other = Constant(other) if not other: return self return Sum(*(self._Children + (-other,))) @@ -288,23 +280,19 @@ class Product(Expression): return isinstance(other, Product) and (self._Children == other._Children) def __mul__(self, other): - if not isinstance(other, Expression): - other = Constant(other) - elif isinstance(other, Product): + if isinstance(other, Product): return Product(*(self._Children + other._Children)) if not other: - return Constant(0) + return 0 if not other-1: return self return Product(*(self._Children + (other,))) def __rmul__(self, other): - if not isinstance(other, Expression): - other = Constant(other) - elif isinstance(other, Product): + if isinstance(other, Product): return Product(*(other._Children + self._Children)) if not other: - return Constant(0) + return 0 if not other-1: return self return Product(*((other,) + self._Children)) @@ -418,17 +406,15 @@ def make_variable(var_or_string): def subscript(expression, index): - if not isinstance(index, Expression): - index = Constant(index) return Subscript(expression, index) -def sum(*components): +def sum(components): components = tuple(c for c in components if c) if len(components) == 0: - return Constant(0) + return 0 elif len(components) == 1: return components[0] else: @@ -438,9 +424,9 @@ def sum(*components): def linear_combination(coefficients, expressions): - return sum(*[coefficient * expression + return sum(coefficient * expression for coefficient, expression in zip(coefficients, expressions) - if coefficient and expression]) + if coefficient and expression) @@ -448,12 +434,12 @@ def linear_combination(coefficients, expressions): def product(*components): for c in components: if not c: - return Constant(0) + return 0 components = tuple(c for c in components if (c-1)) if len(components) == 0: - return Constant(1) + return 1 elif len(components) == 1: return components[0] else: @@ -469,15 +455,17 @@ def polynomial_from_expression(expression): def quotient(numerator, denominator): - if not isinstance(numerator, Expression): - numerator = Constant(numerator) - if not isinstance(denominator, Expression): - denominator = Constant(denominator) - + if not (denominator-1): + return numerator + + import pymbolic.rational as rat + if isinstance(numerator, rat.Rational) and \ + isinstance(denominator, rat.Rational): + return numerator * denominator.reciprocal() + try: - import pymbolic.rational as rat - if isinstance(traits.common_traits(numerator, denominator), - traits.EuclideanRingTraits): + c_traits = traits.common_traits(numerator, denominator) + if isinstance(c_traits, traits.EuclideanRingTraits): return rat.Rational(numerator, denominator) except traits.NoCommonTraitsError: pass diff --git a/src/rational.py b/src/rational.py index 8c42a90..17a4ce0 100644 --- a/src/rational.py +++ b/src/rational.py @@ -65,8 +65,13 @@ class Rational(prm.Expression): gcd_1 = t.gcd(self.Numerator, other.Denominator) gcd_2 = t.gcd(other.Numerator, self.Denominator) - return Rational(self.Numerator/gcd_1 * other.Numerator/gcd_2, - self.Denominator/gcd_2 * other.Denominator/gcd_1) + new_num = self.Numerator/gcd_1 * other.Numerator/gcd_2 + new_denom = self.Denominator/gcd_2 * other.Denominator/gcd_1 + + if not (new_denom-1): + return new_num + + return Rational(new_num, new_denom) except traits.NoCommonTraitsError: return prm.Expression.__mul__(self, other) @@ -90,6 +95,12 @@ class Rational(prm.Expression): def __float__(self): return float(self.Numerator) / flaot(self.Denominator) + def __getinitargs__(self): + return (self.Numerator, self.Denominator) + + def reciprocal(self): + return Rational(self.Denominator, self.Numerator) + def invoke_mapper(self, mapper, *args, **kwargs): return mapper.map_rational(self, *args, **kwargs) diff --git a/src/traits.py b/src/traits.py index b77c972..8c0b59e 100644 --- a/src/traits.py +++ b/src/traits.py @@ -16,8 +16,8 @@ def traits(x): try: return x.traits() except AttributeError: - if isinstance(x, (complex, float)): return FieldTraits - if isinstance(x, int): return IntegerTraits + if isinstance(x, (complex, float)): return FieldTraits() + if isinstance(x, int): return IntegerTraits() raise NoTraitsError @@ -102,3 +102,4 @@ class IntegerTraits(EuclideanRingTraits): return 1 else: raise RuntimeError, "0 does not have a prime factor decomposition" + -- GitLab