DSM's profile - activity

Say you have some formula in the form of an expression:

sage: x,y,a,b = var("x y a b")
sage: f = 10*x+2*y
sage: type(f)
<type 'sage.symbolic.expression.Expression'>

Now I want to replace the 10 by the variable a. What's the easiest way to do this? It's the opposite of the typical substitution, so the functions I'd usually throw at it all fail.

So far I can only up with (1) string techniques, which are bugs waiting to happen, or (2) walking the entire expression tree and constructing a new expression from each operator/operand triplet.

Your error message says `No module named scilab`, not `No module named scipy`, so there was a typo somewhere.

Updated: in modern Sage, alarm raises an AlarmInterrupt, not a KeyboardInterrupt, and it's probably a better idea to use cancel_alarm() instead of alarm(0), which was apparently never offically supported.

Yes, you can use the alarm function and catch the resulting KeyboardInterrupt

sage: alarm(3)
[3 seconds later..]
sage:  
KeyboardInterrupt

Here's an example use:

def long_calculation(sleepfor):
    alarm(3)
    try:
        print 'starting calculation..'
        sleep(sleepfor)
    except KeyboardInterrupt:
        print 'did not complete!'
    # if the computation finished early, though, the alarm is still ticking!
    # so let's turn it off..
    alarm(0) 

sage: long_calculation(6)
starting calculation..
did not complete!
sage: long_calculation(2)
starting calculation..
sage:

But this should simply work, you shouldn't need to play timing games. After interrupting one of the really long calculations, I managed to break an assert:

sage: eq=1/2*((2*e^(2*x) + e^(4*x) + 1)*(e^(2*x) + 1)^(1/2*e^(-x))*2^(1/2*e^x)*e^(1/2*x*e^x)*log(2*e^x/(e^(2*x) + 1)) + (2*e^(2*x) - e^(4*x) - 1)*(e^(2*x) + 1)^(1/2*e^(-x))*2^(1/2*e^x)*e^(1/2*x*e^x))*e^(-1/2*x*e^(-x))/((e^(3*x) + e^x)*(e^(2*x) + 1)^(1/2*e^x)*2^(1/2*e^(-x))) == 1/2*(2*(2*e^(2*x) + e^(4*x) + 1)*e^(x^(-1/2*e^(-x) + 1/2*e^x))*log(2*e^x/(e^(2*x) + 1)) + (2*e^(2*x) - e^(4*x) - 1)*e^(2*x^(-1/2*e^(-x) + 1/2*e^x)) + 2*e^(2*x) - e^(4*x) - 1)/((e^(3*x) + e^x)*e^(2*x^(-1/2*e^(-x) + 1/2*e^x)) + e^(3*x) + e^x)
sage: bool(eq)
get_z.c:37: MPFR assertion failed: (!(((r)->_mpfr_exp) == (((-9223372036854775807L - 1L))+2)) && !(((r)->_mpfr_exp) == (((-9223372036854775807L - 1L))+3)))
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)

so clearly things are getting themselves into a state.

UPDATE: I should also give the usual warning about using "bool(eq)", which is that it doesn't do what you might think it does. bool(eq) == False does not mean that the two sides aren't equal; it could mean only that it couldn't figure out how to prove they were equal. I disagree with this (maxima-based, I think) design decision, and would return True, False, or some Undecided result -- possibly the original expression.

Updated: modern Sage raises an AlarmInterrupt instead of a KeyboardInterrupt, so you can catch it explicitly instead. That's a much cleaner solution than having to dig around inside a KeyboardInterrupt message.

Also, it's probably a better idea to use cancel_alarm() than alarm(0), which was apparently never officially supported.

Sure. While I would have made it a little more explicit, the KeyboardInterrupt generated by alarm does contain information that a typical KeyboardInterrupt doesn't, in the .args and .message attributes. For example, control-C in the following

try:
    sleep(5)
except KeyboardInterrupt as KI:
    print 'args=',KI.args
    print 'message=', KI.message

produces

^Cargs= ()
message=

but adding alarm(2) before executing it produces:

args= ('computation timed out because alarm was set for 2 seconds',)
message= computation timed out because alarm was set for 2 seconds

So one way to take advantage of this would be through something like this:

while True:
    alarm(3)
    try:
        print 'sleeping..'
        sleep(5)
    except KeyboardInterrupt as KI:
        if KI.message.startswith("computation timed out"):
            print 'timed out!'
            continue
        else:
            # not a timeout, so break
            print 'manual exit!'
            alarm(0)
            break
    alarm(0)

which gives

sleeping..
timed out!
sleeping..
timed out!
sleeping..
timed out!
sleeping..
^Cmanual exit!

I tend to wrap this stuff up in decorators, but in this case mine is nested three deep so it's probably overkill for the problem, because it makes it look more complicated than it really is.

No need to use indices; `[f.function(x) for f in F]` should work too.

The Legendre polynomial itself looks fine. To avoid the numerical imprecisions, you could work over QQ instead:

sage: legendre_P(90,t)(0.9)
1.56319477427853e13
sage: legendre_P(90,t)(QQ(0.9)).n()
-0.0874931768791928

It looks like you need 100+ bits of precision in order to evaluate the function in the original order:

sage: legendre_P(90, t)(RealIntervalField(100)("0.9"))
0.?e1
sage: legendre_P(90, t)(RealIntervalField(120)("0.9"))
-0.08749?
sage: legendre_P(90, t)(RealIntervalField(150)("0.9"))
-0.08749317687919?

The polynomial is very oscillatory in that there's a large amount of cancellation going on.

If you're really interested in the value, and only secondarily the polynomial, then you could use either scipy or mpmath, depending on how much precision you'll need:

sage: import scipy.special
sage: scipy.special.eval_legendre(90, 0.9)
-0.08749317687919235
sage: import mpmath
sage: mpmath.legendre(90, 0.9)
mpf('-0.08749317687919235')

@kcrisman: this might be a questiong about the Other Sage..

[spam message]

[spam message removed]

Could you edit your question to include the smallest example you can find which reproduces the problem? Making one will either help you discover the problem yourself or give us information to work with.

That's what you'd get if `sqrt` were actually `math.sqrt` and not Sage's `sqrt`.

Are you sure there are any such pairs to find? Unless I'm missing something, you can't fit $2p^2$ into $(a+1) (b+1) (a+b+2)$ successfully.

None of the Cythonic stuff would work, I don't think, which rules out a lot of the core.

@LuizRobertoMeier: I deleted your insulting comment. When everyone else except you is finding you unclear, you might want to consider the possibility that it's not that everyone except you is an idiot, it's that you've been more ambiguous than you think.

Oh, that's nice.

It's not completely clear to me what you mean by "indexing" a graph edge. If you want to give them labels, you can do that using set_edge_label:

sage: D = DiGraph({ 0: [1,2,3], 1: [0,2], 2: [3], 3: [4], 4: [0,5], 5: [1] })
sage: D.edges()
[(0, 1, None),
 (0, 2, None),
 (0, 3, None),
 (1, 0, None),
 (1, 2, None),
 (2, 3, None),
 (3, 4, None),
 (4, 0, None),
 (4, 5, None),
 (5, 1, None)]
sage: for i, (u, v, l) in enumerate(D.edges()):                                                            
....:     D.set_edge_label(u, v, "edge # {}".format(i))
....:     
sage: D.edges()
[(0, 1, 'edge # 0'),
 (0, 2, 'edge # 1'),
 (0, 3, 'edge # 2'),
 (1, 0, 'edge # 3'),
 (1, 2, 'edge # 4'),
 (2, 3, 'edge # 5'),
 (3, 4, 'edge # 6'),
 (4, 0, 'edge # 7'),
 (4, 5, 'edge # 8'),
 (5, 1, 'edge # 9')]
sage: D.show(edge_labels=True)

produces

IIUC, the numpy library has a routine to do this already:

sage: L = [2,3,1,6,8,-3,9]
sage: numpy.argsort(L)+1
array([6, 3, 1, 2, 4, 5, 7])

Alternatively, you can use a Schwartzian transform (a.k.a. the "decorate-sort-undecorate" idiom). Starting from our list, we can enumerate it, counting from 1 (the "decorate" part)

sage: list(enumerate(L, 1))
[(1, 2), (2, 3), (3, 1), (4, 6), (5, 8), (6, -3), (7, 9)]

Sort these by the second term, the value:

sage: sorted(enumerate(L, 1), key=lambda x: x[1])
[(6, -3), (3, 1), (1, 2), (2, 3), (4, 6), (5, 8), (7, 9)]

Extract the sorted indices:

sage: ii = [pair[0] for pair in sorted(enumerate(L, 1), key=lambda x: x[1])]
sage: ii
[6, 3, 1, 2, 4, 5, 7]

Finally, we can make this into a proper Permutation:

sage: Permutation(ii)
[6, 3, 1, 2, 4, 5, 7]
sage: Permutation(ii).to_cycles()
[(1, 6, 5, 4, 2, 3), (7,)]

and confirm it by applying it to our list:

sage: P = Permutation(ii)
sage: P.action(L)
[-3, 1, 2, 3, 6, 8, 9]