On Mar 30, 5:36 pm, William Stein <wst...@gmail.com> wrote:
> On Mon, Mar 30, 2009 at 8:08 AM, Maurizio <maurizio.gran...@gmail.com> wrote:
>
> > I would really like to not have to annoy you with this stuff, but I
> > really think I'm missing something important (and useful!!)
>
> > The first thing I have to say is: how do I check which is the type of
> > the coefficients (whether they are rationals or something else)? Even
> > when I do multivariate polynomials, I am comfortable with using
> > symbolic variables (probably I should switch, but I still have to
> > learn how to threat polynomials in SAGE).
>
> > One thing that really surprises me is: how comes in this simple
> > following example, expand() is so much slower than factor()? I often
> > do factor(expand(argument)) to be sure that factor is effective, and I
> > was supposing that expand() should have been so straightforward to not
> > require any particular computational effort! But look at this:
>
> > var('x1 x2')
> > var('u')
> > var('s')
> > var('a1 b1 c1 d1')
> > var('a2 b2 c2 d2')
> > var('e1 f1')
> > var('e2 f2')
>
> > A1 = matrix([[a1,b1],[c1,d1]])
> > A2 = matrix([[a2,b2],[c2,d2]])
> > B1 = matrix([[e1],[f1]])
> > B2 = matrix([[e2],[f2]])
>
> > var('D')
> > Am = D*A1 + (1-D)*A2
> > Bm = D*B1 + (1-D)*B2
> > Id=identity_matrix(parent(s),Am.nrows());
> > solm = (Am-Id*s)\(-Bm*u); solm
>
> > [u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) - (b1*D + b2*(1 -
> > D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D - e2*(1 -
> > D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D) - s)*((b1*D + b2*(1 -
> > D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s) + d1*D + d2*(1 - D) -
> > s))]
> > [
> > (u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D - e2*(1 -
> > D))/(a1*D + a2*(1 - D) - s))/((b1*D + b2*(1 - D))*(-c1*D - c2*(1 -
> > D))/(a1*D + a2*(1 - D) - s) + d1*D + d2*(1 - D) - s)]
>
> > sol1 = solm[0];
> > time a = SR(repr(sol1)).factor()
>
> > Time: CPU 0.04 s, Wall: 0.17 s
>
> > time a = SR(repr(sol1)).expand()
>
> > Time: CPU 69.00 s, Wall: 116.28 s
>
> > Do you consider this acceptable?
>
> Of course not. Note that the above step that takes 116.28s on your
> computer takes 1.5 seconds on my laptop, so there is something weird
> about your computer. Using Pynac takes 0.01 seconds and Sympy takes
> 0.02 seconds.
>
> Using current Maxima-based symbolics:
>
> sage: reset()
> sage: var('x1 x2 u s a1 b1 c1 d1 a2 b2 c2 d2 e1 f1 e2 f2 D')
> (x1, x2, u, s, a1, b1, c1, d1, a2, b2, c2, d2, e1, f1, e2, f2, D)
> sage: sol1 = u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) - (b1*D + b2*(1 -
> D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D - e2*(1 -
> D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D) - s)*((b1*D + b2*(1 -
> D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s) + d1*D + d2*(1 - D) -
> s));
> sage:
> sage: time g = sol1.expand()
> CPU times: user 0.40 s, sys: 0.03 s, total: 0.43 s
> Wall time: 1.59 s
>
I was fearing this happened because of the notebook. I can confirm
that even from the command line, the same commands you typed are
giving me the same results: MAXIMA based expand is taking forever...
Maybe by the time I finish writing this post, I'll tell you how much
time.
Here it is...
CPU times: user 57.80 s, sys: 1.67 s, total: 59.47 s
Wall time: 88.15 s
Please, consider that this is a workstation (not even a laptop) and
I'm running SAGE 3.3 (actually, upgraded from 3.2.3) in a virtual
machine with VMWare.
> Using Pynac:
>
> sage: var('x1 x2 u s a1 b1 c1 d1 a2 b2 c2 d2 e1 f1 e2 f2 D', ns=1)
> (x1, x2, u, s, a1, b1, c1, d1, a2, b2, c2, d2, e1, f1, e2, f2, D)
> sage: sol1 = u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) - (b1*D + b2*(1 -
> ....: D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D - e2*(1 -
> ....: D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D) - s)*((b1*D + b2*(1 -
> ....: D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s) + d1*D + d2*(1 - D) -
> ....: s));
> sage:
> sage: time g = sol1.expand()
> CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
> Wall time: 0.00 s
> sage: len(str(g))
> 9533
>
> Sympy is also fast for this sort of thing:
>
> sage: from sympy import var
> sage: var('x1 x2 u s a1 b1 c1 d1 a2 b2 c2 d2 e1 f1 e2 f2 D')
> (x1, x2, u, s, a1, b1, c1, d1, a2, b2, c2, d2, e1, f1, e2, f2, D)
> sage: x1,x2,u,s,a1,b1,c1,d1,a2,b2,c2,d2,e1,f1,e2,f2,D=var('x1 x2 u s
> a1 b1 c1 d1 a2 b2 c2 d2 e1 f1 e2 f2 D')
> sage: sol1 = u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) - (b1*D + b2*(1 -
> D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D - e2*(1 -
> D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D) - s)*((b1*D + b2*(1 -
> D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s) + d1*D + d2*(1 - D) -
> s));
> sage:
> sage: time g=sol1.expand()
> CPU times: user 0.02 s, sys: 0.00 s, total: 0.02 s
> Wall time: 0.02 s
>
Yes, both Pynac and Sympy are fast in doing this on my computer as
well.
> > I think there is some issue here!
> > Please, notice that this is SAGE 3.3... Is this a good reason to
> > upgrade? This is a multiuser environment, so it is slower to upgrade
> > than my laptop.
>
> Upgrading probably won't help. Us switching Sage over to use Pynac
> will help dramatically as you can see above.
>
> > Another question: is possible that whenever I got Wall time much
> > higher than CPU time in factor, this should be recognized as a symptom
> > of having factor done via maxima?
>
> Yes. Anything done purely in Maxima goes entirely to wall time.
>
>
>
> > E.g.:
>
> > time igrFac = SR(repr(igrSol)).factor()
> > Time: CPU 0.51 s, Wall: 2.07 s
>
> > This is a completely different example. Also in this case, the expand
> > took forever to accomplish...
>
> >> Here's an example where the factorization could be done quickly via
> >> Singular, but isn't because we haven't bothered:
>
> >> sage: g = SR(f).subs(x=pi) # replace x by the transcendental pi
> >> sage: time g.factor()
> >> CPU times: user 0.02 s, sys: 0.01 s, total: 0.03 s
> >> Wall time: 1.55 s
>
> > About this point, does substitution by 0 matter as well to make
> > Singular unusable? This would be really bad, because I usually do
> > something like this:
>
> > simpleDict = {"rQ":0, "rCb":0, "rCf":0,"rLf":0}
> > [soluz[0][vCb].subs(simpleDict)
>
> > to simplify my expressions.
>
> I don't even know what you're asking.
>
>
I just wanted some more information about why we haven't bothered to
make that calculation with Singular. It seems that having done the
substitution with pi, prevented evaluating factor() with Singular, and
forced to do this in Maxima. I'm asking whether doing substitution of
some variable to 0 has the same effect. I hope not.
Moreover, I wanted to know if there's any way to recognize whether
Singular has the chance to perform factor() or not, by checking
something like the type of the involved variables, or stuff like that.
>
>
>
> > Last but not least. Do you think is this a reasonable way to work with
> > the solution of a linear system, represented by matrices? The matrices
> > are Aave, Bave, Cave, Dave (canonical dynamic system representation)
>
> > inVec = matrix([[gr],[ib]]); inVec
> > xVec = matrix([[vCb],[vCf],[iLf]]); xVec
> > sys = (Aave - I*w*Id)*xVec + Bave*inVec
> > xSys = [vCb,vCf,iLf]
> > soluz = solve(sys,xSys, solution_dict = True)
> > xSol = matrix([[soluz[0][vCb].subs(simpleDict)],[soluz[0][vCf].subs
> > (simpleDict)],[soluz[0][iLf].subs(simpleDict)]])
> > outVec = Cave*xSol + Dave*inVec
> > igrSol = SR(repr(outVec[0])).subs(simpleDict)
> > vbSol = SR(repr(outVec[1])).subs(simpleDict)
>
> > I fear that I'm confusing SAGE by going back and forth with different
> > packages, so this is taking me hours to do simple operations (am I
> > wrong) like expand(), factor(), collect() and others (especially
> > expand, as you could see!)
>
> > Many, many thanks'
>
> I'm sure doing things like SR(repr(...)) is a very very bad idea. Why
> are you doing that? E.g., why not just
> igrSol = outVec[0].subs(simpleDict)
>
The reason I'm doing this is that otherwise I cannot do factor() or
other functions, because I get
factor(outVec[0])
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/home/notebook/sage_notebook/worksheets/admin/5/code/64.py",
line 7, in <module>
factor(outVec[_sage_const_0 ])
File "/usr/local/sage/local/lib/python2.5/site-packages/
zope.interface-3.3.0-py2.5-linux-i686.egg/", line 1, in <module>
File "/usr/local/sage/local/lib/python2.5/site-packages/sage/rings/
arith.py", line 1800, in factor
raise TypeError, "unable to factor n"
TypeError: unable to factor n
> Or you could try working with a multivariate polynomial ring, which
> *always* expands everything.
>
> sage: R.<x1,x2,u,s,a1,b1,c1,d1, a2, b2, c2, d2, e1, f1, e2, f2, D> = QQ[]
> sage: timeit('u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) - (b1*D +
> b2*(1 - D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D -
> e2*(1 - D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D) - s)*((b1*D +
> b2*(1 - D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s) + d1*D +
> d2*(1 - D) -s))')
> 125 loops, best of 3: 4.44 ms per loop
>
I'll certainly try to work with polynomials, thanks! Is it possible to
make polynomials interact with matrices and linear systems as well
(like solve and other stuff)?
Thanks
Maurizio
> Compare timing:
>
> PYNAC:
> sage: var('x1 x2 u s a1 b1 c1 d1 a2 b2 c2 d2 e1 f1 e2 f2 D',ns=1)
> (x1, x2, u, s, a1, b1, c1, d1, a2, b2, c2, d2, e1, f1, e2, f2, D)
> sage: timeit('expand(u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) -
> (b1*D + b2*(1 - D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 -
> D))*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D)
> - s)*((b1*D + b2*(1 - D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s)
> + d1*D + d2*(1 - D) -s)))')
> 125 loops, best of 3: 7.77 ms per loop
>
> SYMPY:
> sage: from sympy import var
> sage: x1,x2,u,s,a1,b1,c1,d1,a2,b2,c2,d2,e1,f1,e2,f2,D=var('x1 x2 u s
> a1 b1 c1 d1 a2 b2 c2 d2 e1 f1 e2 f2 D')
> sage: timeit('(u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) - (b1*D
> + b2*(1 - D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D
> - e2*(1 - D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D) - s)*((b1*D
> + b2*(1 - D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s) + d1*D +
> d2*(1 - D) -s))).expand()')
> 5 loops, best of 3: 235 ms per loop
>
> MAXIMA-based Symbolics:
> sage: var('x1 x2 u s a1 b1 c1 d1 a2 b2 c2 d2 e1 f1 e2 f2 D',ns=0)(x1,
> x2, u, s, a1, b1, c1, d1, a2, b2, c2, d2, e1, f1, e2, f2, D)sage:
> timeit('expand(u*(-e1*D - e2*(1 - D))/(a1*D + a2*(1 - D) - s) - (b1*D
> + b2*(1 - D))*(u*(-f1*D - f2*(1 - D)) + u*(-c1*D - c2*(1 - D))*(-e1*D
> - e2*(1 - D))/(a1*D + a2*(1 - D) - s))/((a1*D + a2*(1 - D) - s)*((b1*D
> + b2*(1 - D))*(-c1*D - c2*(1 - D))/(a1*D + a2*(1 - D) - s) + d1*D +
> d2*(1 - D) -s)))')
> 5 loops, best of 3: 1.56 s per loop
>
> One thing to notice above is that though sympy appeared nearly as fast
> as Pynac earlier in this email (when using time), that was because of
> lack of timer resolution. In fact, Pynac is 30 times faster in this
> benchmark.
>
> -- William
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to
sage-devel-unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
