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 ? 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 ? 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. 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. 5. I think coercion is important for matrix algorithms. Other than Poly, which other lower-level classes allow coercion ? - 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.
