Hi, On 29 May 2011 12:44, SherjilOzair <[email protected]> wrote:
> Thanks for the pointers. I'm working on what you said. > > Here are a few questions. > > 1. Poly(123) should not give an error, why not treat it as a constant > polynomial? > In this setup 123 means 123*x**0 or 123*x**0*y**0*z**0? Without additional information which is implied by calling 123 a polynomial, 123 is number. Consider a different example. Suppose M and N are matrices. What is M + N? Or M.det()? M and N have to have shape defined, otherwise M and N are just meaningless symbols. > > 2. I used construct_domain on the list of elements of the Matrix, It > returned DMP type, which does *not* allow coercion. What to do if one > of my algorithms involve a square root ? > You have to change the ground domain, either to EX (the simples solution) or to an algebraic domain. > > 3. Even when I tried operating on DMPs using an algorithm which did > not have a square root, an "unsupported operand type(s) for /: 'int' > and 'DMP'" TypeError was returned. > That's fine, because you you have operands of different type. I could be convenient to define DMP() / int, but this is really not necessary. Assuming that 'domain' is an instance of PolynomialRing, then you can do DMP() / domain(int) or even better domain.quo(DMP(), domain(int)), to make division ground type agnostic. > > 4. Should I write different code for the algorithms for each > groundtype ? For example, when using Sympy's type, using Add(*(...)) > adds efficiency, but it doesn't make sense for other types. I can use > sum(...) but that will sacrifice performance slightly. > If you want to write generic algorithms, you can't assume that you work with any particular ground type. The same as you can't assume that integer multiplication is implemented efficiently (mpz vs. int), the same you can't assume that Add() exists. What exists is + operator. It's SymPy's problem that sum(X) is O(len(X)**2) algorithm, but this is just because we use lists instead of dicts for the internal implementation of Add. This will change in future. Also can you say specifically what you want to "Add", because you may be misusing Add for constructing polynomials? > > 5. I think coercion is important for matrix algorithms. Other than > Poly, which other lower-level classes allow coercion ? > Poly is a high-level class. Most types allow for some sort of coercions, so please specify precisely what kind of coercions for what types do you need. Just remember that sympy.polys doesn't depend directly on built-in coersion rules (see DMP() / int discussion), but implements coercions on top of built-in coercions in domains (see Domain.convert() and all from_ and to_ methods). > > - Sherjil Ozair > > On May 28, 11:34 pm, Mateusz Paprocki <[email protected]> wrote: > > Hi, > > > > On 28 May 2011 15:30, SherjilOzair <[email protected]> wrote: > > > > > Hello everyone, > > > I've been successful in writing the symbolic cholesky decomposition > > > for sparse matrices in O(n * c**2) time. This is a reasonable order > > > for sparse systems, but still the performance is not very good. Using > > > python bulitins, It factors a 100 * 100 Matrix with sparsity 0.57 in > > > 961 milli-seconds. Using Sympy's numbers, it takes forever or is pain- > > > stakingly slow for matrices larger than 20 * 20. > > > > In [1] you will find a very simple comparison of Integer, int and mpz. > This > > applies to the rational case as well, just the difference is even bigger. > > > > > > > > > > > > > > > > > > > > > > > > > I understand why we must integrate groundtypes in matrices to make it > > > usable. But I don't know how exactly to do it. > > > > > I see that the Matrix constructor currently employs sympify, so it > > > changes regular ints to Sympy's Integer. I had removed this when I > > > wanted to test for the python builtins in my DOKMatrix implementation. > > > > > Here's an idea that we can build on. Add a kwarg argument in the > > > Matrix constructor called dtype, which could takes values like 'gmpy', > > > 'python', 'sympy', etc. or more detailed, 'int', 'expr', 'poly' etc.. > > > So that, before putting the element in the matrix, we convert it to > > > the required dtype. eg. val = gmpy.mpz(val) > > > > > Is it as simple as this, or am I missing something ? > > > > Following sympy.polys design means that you have to employ static typing > > (all coefficients in a matrix are of the same type, governed by a domain > > that understands properties of the type). Suppose we have a matrix M over > > ZZ, then M[0,0] += 1 is well defined and is fast because it requires only > > one call to domain.convert() (which will exit almost immediately, > depending > > whether ZZ.dtype is Integer, int, mpz or something else). That was > simple, > > but what about M[0,0] += S(1)/2? 1/2 not in ZZ so += may either fail > because > > there is no way to coerce 1/2 to an integer, but it may also figure out a > > domain for 1/2 (QQ), upgrade the domain in M and proceed. In polys both > > scenarios can happen depending whether you use DMP (or any other > low-level > > polynomial representation) or low-level APIs of Poly or high-level APIs > of > > Poly (low-level uses the former and high-level uses the later). The main > > concern in this case is speed (and type checking but it isn't very > strong). > > Deciding whether 1 is in ZZ is fast, but figuring out a domain for 1/2, > > unifying domain of 1/2 with ZZ (a sup domain has to be found, which in > this > > case is simple, but may be highly non-trivial in case of composite or > > algebraic domains) and coercing all elements of M, is slow. > > > > How to figure out a domain for a set of coefficients? Use > > construct_domain(). It will give you the domain and will coerce all > inputs. > > Refer to Poly.__new__ and all Poly.from_* and Poly._from_* to see how > this > > works. sympy.polys should have all tools you will need, so try not to > > reinvent things that are already in SymPy. For example speaking about > those > > "detailed types" 'int', 'expr', 'poly': poly -> what domain and > variables?, > > expr -> what simplification algorithm?, etc. Learn to use what the > library > > provides to you. If there is something missing, e.g. you would like > > construct_domain() to work with nested lists, that can be done, either on > > your own or just ask it. For now it may be a little tedious to use stuff > > from sympy.polys in matrices and at some point I will have share, e.g. > > domains, with other modules. > > > > My suggestion is to start from something simple. You can create a new > matrix > > class that will support the bare minimum of operations to replace Matrix > in > > solve_linear_system(). This new matrix class would support domain > > construction and type conversions using mechanisms from sympy.polys. > Change > > solve_linear_system() to not use simplify() but rely on the ground types > to > > do the job (solving zero equivalence problem). If this works and is fast, > > then you can build on top of this. > > > > > I would like Mateusz especially to comment on this, and also guide me > > > and help me learn about the Polys structure. > > > > [1] > http://mattpap.github.com/masters-thesis/html/src/internals.html#benc... > > > > > > > > > Regards, > > > Sherjil Ozair > > > > > -- > > > 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. > > > > Mateusz > > -- > 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. > > Mateusz -- 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.
