Ask Your Question

DSM's profile - activity

2019-10-18 08:44:31 -0500 received badge  Great Answer (source)
2018-06-06 21:04:22 -0500 received badge  Nice Answer (source)
2018-04-27 14:54:49 -0500 received badge  Good Answer (source)
2017-12-28 02:31:02 -0500 received badge  Great Answer (source)
2017-03-30 10:26:52 -0500 received badge  Guru (source)
2017-03-30 10:26:52 -0500 received badge  Great Answer (source)
2017-02-13 10:06:41 -0500 received badge  Nice Answer (source)
2016-10-07 22:37:17 -0500 received badge  Great Answer (source)
2016-08-06 03:05:58 -0500 received badge  Guru (source)
2016-08-06 03:05:58 -0500 received badge  Great Answer (source)
2016-07-22 17:23:19 -0500 received badge  Famous Question (source)
2016-06-27 18:45:19 -0500 received badge  Good Answer (source)
2016-05-30 17:05:16 -0500 received badge  Nice Answer (source)
2015-09-26 11:51:01 -0500 received badge  Famous Question (source)
2015-02-16 08:00:52 -0500 received badge  Nice Answer (source)
2014-12-11 07:04:13 -0500 received badge  Nice Answer (source)
2014-07-24 07:11:00 -0500 received badge  Good Answer (source)
2014-07-24 07:10:25 -0500 received badge  Nice Answer (source)
2014-06-28 20:14:48 -0500 marked best answer substituting expressions for numbers

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.

2014-01-31 06:26:04 -0500 commented question "No module named scipy"

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

2014-01-27 10:23:42 -0500 edited answer Any way to timeout?

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.

2014-01-27 10:21:47 -0500 edited answer Distinguish between alarm() and manual interrupt

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.

2014-01-16 10:03:38 -0500 commented answer Function generator

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

2014-01-16 07:11:24 -0500 answered a question Bug in Legendre polynomial.

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')
2013-11-25 14:14:44 -0500 received badge  Notable Question (source)
2013-10-04 10:17:50 -0500 commented question Changing country name

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

2013-09-17 01:07:42 -0500 received badge  Nice Answer (source)
2013-09-14 19:41:37 -0500 received badge  Good Answer (source)
2013-08-28 18:10:51 -0500 edited answer [spam message]

[spam message]

2013-08-28 18:10:33 -0500 edited question [spam message]

[spam message removed]

2013-08-16 04:10:53 -0500 commented question RAM memory becomes full while running script

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.

2013-08-15 07:55:35 -0500 commented question Computing derivatives

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

2013-08-05 09:01:06 -0500 received badge  Nice Answer (source)
2013-07-08 05:43:47 -0500 received badge  Good Answer (source)
2013-07-08 02:09:49 -0500 received badge  Good Answer (source)
2013-05-23 07:56:36 -0500 commented question How to find instances where $d(a,b) = p^2$ for $p$ a prime

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.

2013-05-04 06:24:49 -0500 received badge  Enthusiast
2013-04-18 10:32:28 -0500 received badge  Great Question (source)
2013-04-02 11:34:21 -0500 commented question Sage + IronPython

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

2013-03-29 13:00:08 -0500 received badge  Good Answer (source)
2013-03-29 03:43:14 -0500 commented question Table of primes

@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.

2013-03-28 10:23:31 -0500 commented answer Finding the permutation that sorts a list

Oh, that's nice.

2013-03-28 07:50:44 -0500 answered a question code for bipartite graphs

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

image description

2013-03-28 06:07:51 -0500 received badge  Great Answer (source)
2013-03-28 06:07:51 -0500 received badge  Guru (source)
2013-03-27 22:31:22 -0500 received badge  Good Answer (source)
2013-03-27 16:38:34 -0500 edited answer Finding the permutation that sorts a list

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]
2013-03-27 15:43:05 -0500 received badge  Nice Answer (source)