OK, that makes sense. So let's go back to this idea of having Add, Mul, and
Pow classes.
How do Add and Mul figure out how to intelligently combine expressions, such
as Add(x, x) -> Mul(2, x)? For example, here's a sample REPL session:
>>> from sympy import *
>>> x = Symbol("x")
>>> Add(x, sin(x), x, sin(x), 3)
3 + 2*x + 2*sin(x)
>>> Add(x, sin(x), x, sin(x), 3).is_Add
True
>>> Add(x, sin(x), x, sin(x), 3).is_Mul
False
>>> Add(x, sin(x), x, sin(x), 3).args
(3, 2*x, 2*sin(x))
>>> Add(x, sin(x), x, sin(x), 3).args[1]
2*x
>>> Add(x, sin(x), x, sin(x), 3).args[1].is_Mul
True
Somehow Add must be figuring out that sin(x) and sin(x) are the same. Maybe
that's happening in flatten? It's kind of magical.
On Fri, Jun 24, 2011 at 9:16 PM, Aaron Meurer <[email protected]> wrote:
> fdiff makes it easier to define the derivative for functions. If
> f'(x) = g(x), then f.fdiff will return g. Then
> Function._eval_derivative will automatically use that to apply the
> chain rule when differentiating an expression containing f. This
> avoids code duplication of the chain rule across every single
> function.
>
> Aaron Meurer
>
> On Fri, Jun 24, 2011 at 9:56 PM, Jeff Pickhardt <[email protected]>
> wrote:
> > Yeah, I see the advantage in that. That's really cool.
> >
> > Here's another one. I see the Derivative class. But cos already has
> fdiff
> > and an _eval_derivative from class inheritance. (Not sure what the
> > difference between the two methods is, but anyway) So why is a
> Derivative
> > class needed? Couldn't you just make it a function, something like:
> >
> > @sympify_it
> > def derivative(f, with_respect_to):
> > return f.fdiff(with_respect_to)
> >
> > The Derivative class has extra methods like _eval_nseries, but wouldn't
> > those just get subsumed with whatever a derivative function returns? For
> > example, if you call fdiff on cos, it returns -sin, which should have a
> > _eval_nseries method anyhow. What's the point of having a separate
> > Derivative class?
> >
> >
> > --
> > Jeff Pickhardt
> > (650) 530-0036
> > [email protected]
> >
> >
> > On Thu, Jun 23, 2011 at 11:55 PM, Aaron Meurer <[email protected]>
> wrote:
> >>
> >> On Fri, Jun 24, 2011 at 12:33 AM, Jeff Pickhardt <[email protected]>
> >> wrote:
> >> > Thanks guys, this is great stuff. You might want to memorialize this
> on
> >> > a
> >> > wiki or something for newcomers to SymPy.
> >>
> >> That's a good idea.
> >>
> >> >
> >> > So a little bit about my background with Python: I'm a pretty
> >> > experienced
> >> > Python programming, I understand operator overloading, I see your
> >> > @_sympifyit decorators all over the place (I'm guessing they're to
> >> > convert 2
> >> > to Integer(2), as Aaron explained). In fact, I started writing my own
> >> > Python symbolics math library years ago.
> >>
> >> OK. I just wanted to be sure to mention it, since to people who don't
> >> know about operator overloading, it can be the most mysterious part of
> >> the whole thing (how does it "know" to convert x + x into 2*x?)
> >>
> >> >
> >> > One follow-up question is:
> >> >
> >> > Is there an advantage to having your own classes for Pow, Add, Mul,
> and
> >> > other operators? Why didn't you just absorb those in classes for
> >> > expressions, functions, etc, and have an expression + expression
> return
> >> > an
> >> > expression?
> >>
> >> It makes for much better object oriented programming. You can see all
> >> the methods of Add, Mul, and Pow specifically handle the behavior on
> >> that type.
> >>
> >> Let's use differentiation as an example. diff(expr, x) is implemented
> >> as expr._eval_derivative(x) (technically, this is not 100% true
> >> because there's also fdiff, but let's suppose that it's just done this
> >> way for simplicity). If you do it this way, it's a very simple
> >> recursive algorithm. You define your base cases, and define
> >>
> >> Add
> >>
> >> def _eval_derivative(self, x):
> >> return Add(*[i._eval_derivative(x) for i in self.args])
> >>
> >> Mul
> >>
> >> def _eval_derivative(self, x):
> >> return Add(*(Mul(*(self.args[:i] +
> >> [self.args[i]._eval_derivative(x)] + self.args[i+1:]))))
> >>
> >> Pow
> >>
> >> def _eval_derivative(self, x):
> >> # I hope I have the rule correct here
> >> return self.exp*self.base._eval_derivative(x)*self.base**(self.exp
> >> - 1) + log(self.base)*self.exp._eval_derivative(x)*self.base**self.exp
> >>
> >> (note, I didn't copy this from the actual SymPy code, so it might be a
> >> little different, but it will be similar).
> >>
> >> This way, the code is all very simple, because each class only needs
> >> to know how to deal with itself.
> >>
> >> This also makes it much more modular. If you want to extend SymPy,
> >> you just create your own class. You can make that class work with
> >> diff() just by defining ._eval_derivative() on that class. With your
> >> method, if someone wants to extend the code, they have to extend the
> >> Expression class (which will be huge, btw).
> >>
> >> Your method goes half way, because it has a separate "Add" class,
> >> PolyTerm, but you are still keeping things like x and x**2 and x*y as
> >> the same class, when the first should be a Symbol, the second a Pow,
> >> and the third a Mul.
> >>
> >> By the way, if you are interested in other ways to implement symbolics
> >> in Python, you might look at the sympycore project, which is a
> >> research project that tries to implement symbolics in the most
> >> efficient way possible (and is often faster than SymPy as a result).
> >> See http://code.google.com/p/sympycore/.
> >>
> >> Aaron Meurer
> >>
> >> >
> >> > Here's an example. Specifically, x^2, instead of being Pow(x, 2),
> would
> >> > be
> >> > something like:
> >> > MonoTerm (self.units = {x: 1}).__pow__(2) becomes MonoTerm, self.units
> =
> >> > {x:
> >> > 2}.
> >> >
> >> > For example, digging up my old code, when I was making a symbolics
> math
> >> > library, I started doing it this way:
> >> >
> >> > class MonoTerm:
> >> > """A MonoTerm is an instance of terms multiplied together, but
> >> > without
> >> > any sums or differences. As an example, 2 could be a MonoTerm, as
> could
> >> > 3
> >> > z, but 2 x + 3 y must be a PolyTerm. MonoTerms can only be operated
> on
> >> > by
> >> > MonoTerms with the same dimensional units."""
> >> > def __init__(self, string=''):
> >> > # instance variables
> >> > self.number=0
> >> > self.units={}
> >> > # initialize to the string, with the become method
> >> > self.become(string)
> >> >
> >> > def become(self, string=''):
> >> > """Handles the parsing of the MonoTerm into number and
> units."""
> >> > if string=='':
> >> > # nothing should happen; this is the "empty term"
> >> > return
> >> > if m.search('^\s*('+patterns.number+')', string):
> >> > self.number = float(m[0])
> >> > else:
> >> > self.number = 1.0 # 'kg' == '1.0 kg'
> >> > unitSearch =
> >> >
> >> >
> re.findall('('+patterns.unit+')\s*(?:(?:\^|\*\*)\s*('+patterns.number+'))?',
> >> > string)
> >> > for unit, exponent in unitSearch:
> >> > if unit == '':
> >> > continue
> >> > if exponent == '':
> >> > self.units[unit]=1.0
> >> > else:
> >> > self.units[unit]=float(exponent)
> >> >
> >> > def __add__(self, other):
> >> > """Adds two MonoTerms together, provided they have the same
> >> > units."""
> >> > if self.units != other.units:
> >> > # can only subtract MonoTerms with MonoTerms
> >> > raise Exception('Cannot add two MonoTerms with
> incompatible
> >> > units. Use PolyTerms instead.')
> >> > try:
> >> > newTerm = MonoTerm()
> >> > newTerm.number = self.number + other.number
> >> > newTerm.units = self.units
> >> > newTerm.prune()
> >> > return newTerm
> >> > except:
> >> > return other+self
> >> > (etc)
> >> >
> >> > class PolyTerm:
> >> > """A PolyTerm is a term made up of the sum of two or more
> >> > MonoTerms."""
> >> > def __init__(self, string=''):
> >> > # instance variables
> >> > self.terms={}
> >> > self.become(string)
> >> >
> >> > def become(self, string=''):
> >> > # assumes the string is well-formed various terms split by +
> >> > # Actually, for now, assume the string is really a MonoTerm...
> >> > this
> >> > will be updated later.
> >> > newTerm = MonoTerm(string)
> >> > self.terms[newTerm.getLabel()]=newTerm
> >> >
> >> > def __add__(self, other):
> >> > """Adds two PolyTerms together."""
> >> > newTerm = copy.deepcopy(self)
> >> > for term in other:
> >> > if term.getLabel() in newTerm.terms:
> >> > newTerm.terms[term.getLabel()] =
> >> > newTerm.terms[term.getLabel()] + term
> >> > else:
> >> > newTerm.terms[term.getLabel()] = copy.deepcopy(term)
> >> > newTerm.prune()
> >> > return newTerm
> >> > (etc)
> >> >
> >> > term1 = PolyTerm('1 x y')
> >> > term2 = PolyTerm('1 y z')
> >> > term3 = PolyTerm('2 x y^2 z')
> >> > term4 = PolyTerm('1 y^2')
> >> >
> >> > print ((term1+term2)**2 - term3) / term4 # ((xy+yz)^2 - 2xy^2z)/y^2 ==
> >> > x^2 +
> >> > z^2
> >> >
> >> >
> >> > --
> >> > Jeff Pickhardt
> >> > (650) 530-0036
> >> > [email protected]
> >> >
> >> >
> >> > On Thu, Jun 23, 2011 at 10:33 PM, Aaron Meurer <[email protected]>
> >> > wrote:
> >> >>
> >> >> Well, you picked a particularly tricky example because it uses both
> >> >> the main core and the polys.
> >> >>
> >> >> So let's start with my_expresion (the core part). You already know
> >> >> that x is defined to be a Symbol object (because you did this
> >> >> yourself). You should also know that sin and cos are Function
> >> >> objects, which you can think of as container objects that hold
> >> >> arguments and which also have various mathematical relations defined
> >> >> on them.
> >> >>
> >> >> As you may know, Python let's you override the behavior of the
> >> >> built-in operations *, +, /, -, etc. on your own objects. So all
> >> >> objects have method __mul__, __add__, etc. defined on them. When you
> >> >> call x + y, this reduces to x.__add__(y). In SymPy, a.__add__(b) is
> >> >> converted to Add(a, b). The same is true for __mul__ and Mul. So in
> >> >> sin(x)**2 + 2*sin(x) + 1, sin(x) creates a sin object (when it is
> >> >> called). This reduces to
> >> >>
> >> >> sin(x).__pow__(2) + sin(x).__rmul__(2) + 1
> >> >> Pow(x, 2).__add__(Mul(2, sin(x))).__add__(1)
> >> >> Add(Pow(x, 2), Mul(2, sin(x)), 1)
> >> >>
> >> >> which is what you will get if you call srepr(my_expression). Note
> >> >> that 2*sin(x) actually calls __rmul__. That's because 2 (type int)
> >> >> doesn't know how to multiply sin(x) (type sin), so sin(x)'s __rmul__
> >> >> method is called. This brings up an important point. All Python
> >> >> types are converted in this process to SymPy types, through the
> >> >> sympify() function. So, for example, sin(x).__pow__(2) reduces to
> >> >> sin(x).__pow__(sympify(2)), which results in
> >> >> sin(x).__pow__(Integer(2)).
> >> >>
> >> >> Now, Mul, Pow, and Add's __new__ methods contain logic to do
> automatic
> >> >> simplifications like x + 2*x => 3*x or x*x => x**2. In this case, no
> >> >> such simplifications needed to be applied.
> >> >>
> >> >> All the args for any SymPy expression are stored in expr.args. So
> you
> >> >> would have
> >> >>
> >> >> >>> my_expression.args
> >> >> (1, sin(x)**2, 2*sin(x))
> >> >> >>> my_expression.args[1].args
> >> >> (sin(x), 2)
> >> >>
> >> >> If you want to see the code for all of this, you can look at the
> files
> >> >> in sympy/core. The __op__ stuff is mostly in basic.py and expr.py.
> >> >> The flattening routines are in add.py, pow.py, and mul.py. The code
> >> >> for functions that sin is built off of is in function.py.
> >> >>
> >> >> Now, to the factoring part. my_expression.factor() is a shortcut to
> >> >> factor(my_expression). Because factoring is a polynomial algorithm,
> >> >> the expression has to be converted to a Poly first. Poly is able to
> >> >> represent polynomials with arbitrary symbolic "generators". In this
> >> >> case, it determines that it should use sin(x) as a generator, so it
> >> >> creates Poly(sin(x)**2 + 2*sin(x) + 1, sin(x)), which you can think
> of
> >> >> as a wrapper around the polynomial y**2 + 2*y + 1, where y is set to
> >> >> be sin(x) (for the purposes of Poly, it does not matter what the
> >> >> generators are, other than that coefficients cannot contain symbols
> >> >> from them, so you can think of it in this way).
> >> >>
> >> >> Then, it calls the factorization algorithm on Poly. If you are
> >> >> interested in how this works, I suggest you read the code and the
> >> >> papers referenced there. In this case, it is able to factor the
> >> >> polynomial using the squarefree factorization algorithm, which is
> >> >> actually not too difficult to understand. In the general case, it
> >> >> uses a complicated multivariate factorization algorithm that factors
> >> >> any multivariate polynomial into irreducibles.
> >> >>
> >> >> Anyway, Poly.factorlist returns something like [(Poly(sin(x) + 1),
> >> >> 2)]. factor() converts this into a normal SymPy expression (also
> >> >> called a Basic expression or Expr expression) by passing it to Mul
> and
> >> >> Pow (something like Mul(*[Pow(b, e) for b, e in expr]) would do it, I
> >> >> think).
> >> >>
> >> >> All the Poly stuff lives in sympy/polys. If you are interested, I
> can
> >> >> explain a little how they work (internal representation, etc.). The
> >> >> code for the Poly class lives in polytools.py, though the actual
> >> >> factorization algorithm lives in sqftools.py and factortools.py.
> >> >>
> >> >> And one thing that I didn't mention (and maybe you didn't even think
> >> >> of) is the printing. You do not use pretty printing, so the printing
> >> >> is rather simple (just recursively print the objects using the proper
> >> >> operators). If you are interested, you can look at the code in
> >> >> sympy/printing.
> >> >>
> >> >> Let me know if this made sense, and if there are bits that you still
> >> >> would like to know about. Also, remember that SymPy is written in a
> >> >> fairly modular way, so it's completely unnecessary to know how a
> >> >> module works unless you want to work on that module specifically
> >> >> (e.g., you don't need to how the core works to write some
> >> >> simplification algorithm, like in simplify.py).
> >> >>
> >> >> Aaron Meurer
> >> >>
> >> >> On Thu, Jun 23, 2011 at 10:42 PM, Jeff Pickhardt <
> [email protected]>
> >> >> wrote:
> >> >> > Hey guys,
> >> >> >
> >> >> > I'm reading through the SymPy code and, well, it's somewhat
> >> >> > overwhelming
> >> >> > if
> >> >> > you're new to the project because there's so much going on.
> (That's
> >> >> > a
> >> >> > good
> >> >> > thing too - it means its robust!)
> >> >> >
> >> >> > Can someone help explain how this works?
> >> >> >
> >> >> >>>> from sympy import *
> >> >> >>>> x = Symbol("x")
> >> >> >>>> my_expression = sin(x)**2 + 2*sin(x) + 1
> >> >> >>>> my_expression.factor()
> >> >> > (1 + sin(x))**2
> >> >> >>>>
> >> >> >
> >> >> > For instance, what data structures happen when I create
> >> >> > my_expression,
> >> >> > what
> >> >> > happens when I factor it, etc. A high-level walk through would
> >> >> > help. I
> >> >> > see
> >> >> > there's stuff going on at polytools.py, and I think
> _symbolic_factor
> >> >> > gets
> >> >> > called. It's just confusing to keep everything in my head when I
> >> >> > don't
> >> >> > yet
> >> >> > have a high level understanding of how sympy expressions and what
> not
> >> >> > actually work.
> >> >> >
> >> >> > Also, the new quantum mechanics stuff looks really cool. Wish I
> had
> >> >> > had
> >> >> > that a few years ago!
> >> >> >
> >> >> > Thanks,
> >> >> > Jeff
> >> >> >
> >> >> > --
> >> >> > Jeff Pickhardt
> >> >> > (650) 530-0036
> >> >> > [email protected]
> >> >> >
> >> >> > --
> >> >> > You received this message because you are subscribed to the Google
> >> >> > Groups
> >> >> > "sympy" group.
> >> >> > To post to this group, send email to [email protected].
> >> >> > To unsubscribe from this group, send email to
> >> >> > [email protected].
> >> >> > For more options, visit this group at
> >> >> > http://groups.google.com/group/sympy?hl=en.
> >> >> >
> >> >>
> >> >> --
> >> >> You received this message because you are subscribed to the Google
> >> >> Groups
> >> >> "sympy" group.
> >> >> To post to this group, send email to [email protected].
> >> >> To unsubscribe from this group, send email to
> >> >> [email protected].
> >> >> For more options, visit this group at
> >> >> http://groups.google.com/group/sympy?hl=en.
> >> >>
> >> >
> >> > --
> >> > You received this message because you are subscribed to the Google
> >> > Groups
> >> > "sympy" group.
> >> > To post to this group, send email to [email protected].
> >> > To unsubscribe from this group, send email to
> >> > [email protected].
> >> > For more options, visit this group at
> >> > http://groups.google.com/group/sympy?hl=en.
> >> >
> >>
> >> --
> >> You received this message because you are subscribed to the Google
> Groups
> >> "sympy" group.
> >> To post to this group, send email to [email protected].
> >> To unsubscribe from this group, send email to
> >> [email protected].
> >> For more options, visit this group at
> >> http://groups.google.com/group/sympy?hl=en.
> >>
> >
> > --
> > You received this message because you are subscribed to the Google Groups
> > "sympy" group.
> > To post to this group, send email to [email protected].
> > To unsubscribe from this group, send email to
> > [email protected].
> > For more options, visit this group at
> > http://groups.google.com/group/sympy?hl=en.
> >
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To post to this group, send email to [email protected].
> To unsubscribe from this group, send email to
> [email protected].
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.
>
>
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
