Re: [sage-support] Reduce all coefficients of multivariate polynomial to 1

[D. S. McNeil]

In Python generally, "X = Y" doesn't modify the *object* which is named X,
it just rebinds the name, saying that "X" now refers to the object given by
the expression Y.  So in general,

for something in someloop:
 something = something_else

isn't going to do much.  (I want to say "nothing" but I guess there could
be side effects of computing something_else.)


You can arrive at your goal in lots of ways.  You could sum the component
monomials, which is short but less general:

sage: R. = QQbar[]
sage: P = x1^2 + 2*x1*x2 + x2^2 + 2*x1*x3 + 2*x2*x3 + x3^2 + 2*x1 + 2*x2 +
2*x3 + 1
sage: sum(P.monomials())
x1^2 + x1*x2 + x2^2 + x1*x3 + x2*x3 + x3^2 + x1 + x2 + x3 + 1


Another way to get what you want would be to use map_coefficients, which
applies the function you pass it to all the coefficients, which is clear
but doesn't give you access to the corresponding term:

sage: P = x1^2 + 2*x1*x2 + x2^2 + 2*x1*x3 + 2*x2*x3 + x3^2 + 2*x1 + 2*x2 +
2*x3 + 1
sage: P.map_coefficients(lambda x: 1)
x1^2 + x1*x2 + x2^2 + x1*x3 + x2*x3 + x3^2 + x1 + x2 + x3 + 1


Or in the most general way, you could loop over the coefficients and their
monomials and build a new polynomial from each term:

sage: list(P)
[(1, x1^2),
 (2, x1*x2),
 (1, x2^2),
 (2, x1*x3),
 (2, x2*x3),
 (1, x3^2),
 (2, x1),
 (2, x2),
 (2, x3),
 (1, 1)]
sage: sum(mon for coeff, mon in P)
x1^2 + x1*x2 + x2^2 + x1*x3 + x2*x3 + x3^2 + x1 + x2 + x3 + 1


Doug

-- 
You received this message because you are subscribed to the Google Groups 
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at https://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.

Re: [sage-support] Weird comprehension behavior.

[D. S. McNeil]

Not scoping, typing.

When you do i+1, this is preparsed into

sage: preparse('[len(top_points(i+1)) for i in range(2,10)]')
'[len(top_points(i+Integer(1))) for i in range(Integer(2),Integer(10))]'

and the addition means that top_points is passed not a Python int, but a
Sage Integer, and your divisions aren't Python 2 truncating division.

To avoid this, you could use srange:

sage: [len(top_points(i)) for i in range(3,11)]
[5, 5, 5, 5, 5, 5, 5, 5]
sage: [len(top_points(i)) for i in srange(3,11)]
[30, 55, 91, 140, 204, 285, 385, 506]

or the `..` notation:

sage: [len(top_points(i)) for i in (3..10)]
[30, 55, 91, 140, 204, 285, 385, 506]

or anything which makes sure that you're not dividing an int by an int.


Doug

-- 
You received this message because you are subscribed to the Google Groups 
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at https://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.

Re: [sage-support] Plotting a q analogue function as a challenge?

[D. S. McNeil]

While qgamma isn't a native function, there's a qgamma implementation in
mpmath, one of the libraries included in Sage, so:

from mpmath import qgamma
plot(lambda x: qgamma(4,x), (x, 2, 10))

should give you a plot of gamma_(q=4).


Doug

-- 
You received this message because you are subscribed to the Google Groups 
sage-support group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.

Re: [sage-support] Random sum of functions

[D. S. McNeil]

I don't think you need to make an explicit class here.  You can build a
function from within another function, and return that:

sage: def f(n, z):
: return z**n
:
sage: def maker(tup):
: def g(z):
: return sum(abs(f(a_i,z))**2 for a_i in tup)
: return g
:

After which

sage: f(1, 3+2*I)
2*I + 3
sage: maker((1,5,3))
function __main__.g
sage: g = maker((1,5,3))
sage: g(5)
9781275
sage: abs(f(1,5))**2 + abs(f(5,5))**2 + abs(f(3,5))**2
9781275

I used Python functions here because that's what you started with.  A
similar approach would have worked if you were making a Sage function
instead, namely write a function which accepts a tuple and returns a new
function.


Doug


On Wed, Jul 17, 2013 at 5:07 AM, Laurent Decreusefond 
laurent.decreusef...@gmail.com wrote:

 Hi everyone,

 say I have a function of both an integer n and a complex z

 def f(n,z):
  return z**n

 For any tuple of integer (a_1,a_2, ..., a_k), (actually, k and a_i are
 random), I want to form the function

 z - sum_{i=1}^k |f(a_i, z)|^2

 The result must still be a function. I guess it is something which can be
 done by redefining a new class but I'm a totally newbie at object
 programming.
 Can you help me ?

  --
 You received this message because you are subscribed to the Google Groups
 sage-support group.
 To unsubscribe from this group and stop receiving emails from it, send an
 email to sage-support+unsubscr...@googlegroups.com.
 To post to this group, send email to sage-support@googlegroups.com.
 Visit this group at http://groups.google.com/group/sage-support.
 For more options, visit https://groups.google.com/groups/opt_out.




-- 
You received this message because you are subscribed to the Google Groups 
sage-support group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/groups/opt_out.

Re: [sage-support] Evaluating a symbolic expression

[D. S. McNeil]

 I mean, I can of course do f(P1[0],P1[1],P1[2]),
 but this is highly non-elegant.

You can use an asterisk:

sage: f(x,y,z) = x+10*y+100*z
sage: P1 = [2,3,4]
sage: f(*P1)
432

Here, * behaves kind of like a dereferencing operator: f(*(x,y,z)) ==
f(x,y,z).  See this StackOverflow question for a longer explanation:

http://stackoverflow.com/questions/2921847/python-once-and-for-all-what-does-the-star-operator-mean-in-python


Doug

-- 
You received this message because you are subscribed to the Google Groups 
sage-support group.
To post to this group, send email to sage-support@googlegroups.com.
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support?hl=en.

Re: [sage-support] SEGV in mpmath.findroot

[D. S. McNeil]

Hmm.  My backtrace showed:

/home/mcneil/sagedev/sage-5.4.beta0/local/lib/libcsage.so(print_backtrace+0x3b)[0xb6bd7c49]
/home/mcneil/sagedev/sage-5.4.beta0/local/lib/libcsage.so(sigdie+0x17)[0xb6bd7c89]
/home/mcneil/sagedev/sage-5.4.beta0/local/lib/libcsage.so(sage_signal_handler+0x212)[0xb6bd778c]
[0xb772c400]
/home/mcneil/sagedev/sage-5.4.beta0/local/lib/libgmp.so.7(__gmpn_rshift+0x19c)[0xb6b897dc]

and my standard Python version of mpmath, which I don't think is
gmp-aware, produces

/usr/local/lib/python2.7/dist-packages/mpmath/libmp/libelefun.pyc in
exp_fixed(x, prec, ln2)
   1410 return v  n
   1411 else:
- 1412 return v  (-n)
   1413
   1414

OverflowError: long int too large to convert to int


Doug

-- 
You received this message because you are subscribed to the Google Groups 
sage-support group.
To post to this group, send email to sage-support@googlegroups.com.
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support?hl=en.

Re: [sage-support] randint - Maybe one bug !

[D. S. McNeil]

On Wed, Sep 19, 2012 at 12:42 PM, Christophe BAL projet...@gmail.com wrote:
 What I think very confusing is that 1/4 is the Sage division and not the
 Python standard one, so why it would be different for randint ?

It's not Sage division vs. Python division, it's Sage Integers vs. Python ints.

At the Sage console, input lines are preparsed to wrap numbers like 1
to 4 in Integer calls:

sage: 1/4
1/4
sage: preparse(1/4)
'Integer(1)/Integer(4)'

In fact, you could even turn the preparser off, and recover Python's behaviour:

sage: preparser(False)
sage: 1/4
0


Doug

-- 
You received this message because you are subscribed to the Google Groups 
sage-support group.
To post to this group, send email to sage-support@googlegroups.com.
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support?hl=en.

Re: [sage-support] Re: Using find_fit() on lambert_w()

[D. S. McNeil]

Somewhere along the path the second argument to scipy.special.lambertw
(0) is being converted to a float:

 import scipy.special
 scipy.special.lambertw
built-in function lambertw
 scipy.special.lambertw(1, 0)
(0.56714329040978384+0j)
 scipy.special.lambertw(float(1), 0)
(0.56714329040978384+0j)
 scipy.special.lambertw(float(1), float(0))
Traceback (most recent call last):
  File stdin, line 1, in module
  File lambertw.pyx, line 340, in scipy.special.lambertw.lambertw
(scipy/special/lambertw.c:1378)
TypeError: ufunc '' not supported for the input types, and the inputs
could not be safely coerced to any supported types according to the
casting rule 'safe'

If you (1) bypass this by using a Python function and making a number
explicitly, (2) use something like the abs of LambertW, and (3) change
your introduced error to be smaller  -- multiplying by
normalvariate(0, 1) is pretty hard to recover from! -- then it should
work:

##
def g(t,a,b,c):
return abs(a*lambert_w(b*t-c).n())

a = 1.2
b = 0.5
c = 0.1

data = [(i, g(i, a, b, c) + normalvariate(0, 0.05)) for i in
xsrange(0, 4*pi, 0.2)]
var('a, b, c, x')
model = g
fit = find_fit(data, model, parameters=[a,b,c], variables=[x])
##

produces

sage: fit
[a == 1.166896678207367, b == 0.5327270673942398, c == 0.1220748294541312]

which isn't bad.


Doug

-- 
You received this message because you are subscribed to the Google Groups 
sage-support group.
To post to this group, send email to sage-support@googlegroups.com.
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support?hl=en.

Re: [sage-support] Out of memory error when dividing vector by scalar

[D. S. McNeil]

After some digging -- and a fortuitous control-C at the right moment
-- it looks like it's trying to construct the basis for the ambient
free module.  This will be a list of 10^4 vectors, each 10^4 elements
long, so it's not surprising it takes a lot of memory.. adding a print
statement in the loop in FreeModule_ambient.basis seems to confirm
that's the issue.

Unfortunately suggesting how to get around this requires a deeper
understanding of Sage internals than I have. :-/


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] slow symbolic computation

[D. S. McNeil]

On Fri, Jun 22, 2012 at 7:19 AM, David Harvey d.har...@unsw.edu.au wrote:
 --
 | Sage Version 5.0, Release Date: 2012-05-14                         |
 | Type notebook() for the GUI, and license() for information.        |
 --
 sage: time floor(log(4)/log(2))
 2
 Time: CPU 1.07 s, Wall: 1.07 s
 sage: time floor(log(4)/log(2))
 2
 Time: CPU 1.07 s, Wall: 1.07 s

 This seems very slow to me. Am I doing something wrong?

No, this is a weird quirk of the way things are implemented at the
moment which means that integer values are actually the worst case..
See the trac ticket

http://trac.sagemath.org/sage_trac/ticket/12121

or my answer here

http://ask.sagemath.org/question/964/lenlist-ceillog42-bugs

for a more detailed explanation.  Should probably fix this at some point..


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Divisor Riemann Roch

[D. S. McNeil]

 AttributeError: 'ProjectiveCurve_finite_field' object has no attribute
 'riemann_roch_basis'

 why?? help please!

If I understand correctly, as the documentation says:

Currently this only works over prime field and divisors supported
on rational points.

Your F is GF(4), which isn't a prime field, so C doesn't offer the
riemann_roch_basis method.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] gauss_sum() core dump

[D. S. McNeil]

On Fri, Jun 1, 2012 at 10:59 AM, Jeroen Demeyer jdeme...@cage.ugent.be wrote:

 Looks like an infinite recursion, leading to stack exhaustion, leading
 to a SIGSEGV:

I agree with the last two but not the first.  I think it's just a
really long LazyBinop chain, which is why I was having trouble
debugging the problematic part of polynomial_element.pyx:

elif (d  4 or d  5) and self._compiled is None:
while i = 0:
result = result * a + self[i]
i -= 1

Even putting in print statements triggered unlazification, or
whatever the word is (evaluation?), which blew the stack.  For
example:

-
from sage.rings.real_lazy import LazyBinop

def loop(n):
x = RealLazyField()(1)
for i in xrange(n):
x = LazyBinop(RLF, x, x, operator.mul)
-

produces

sage: loop(10)
sage: loop(100)
sage: loop(1000)
sage: loop(1)
sage: loop(10)
sage: loop(100)
/home/mcneil/sagedev/sage-5.1.beta1/spkg/bin/sage: line 326:  8729
Segmentation fault  (core dumped) sage-ipython $@ -i


Doug

PS:

Side issue #1: evaluating LazyBinops  is *really* slow as soon as
there get to be more than a few terms.

Side issue # 2: it looks weird to me that we loop over every item in
the polynomial list rather than over every coefficient in the elif
branch I quoted above.  Having lots of zeroes is a pretty common use
case, I'd've thought.

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] gauss_sum() core dump

[D. S. McNeil]

I can push this back at least to
number_field_morphisms.create_embedding_from_approx: the segfault
seems to happen during the evaluation of the defining polynomial of
the cyclotomic field.

For my 5.0,

sage: CyclotomicField(160400)
[...]
/Applications/sage/spkg/bin/sage: line 312: 18832 Segmentation
fault  sage-ipython $@ -i

is enough to crash it.




Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] using a base other then base 10.

[D. S. McNeil]

The following:

sage: z = pi.n(100)
sage: z.str(base=4)
'3.0210033310202011220300203103010301212022023200'

should get you started.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Sage 5.0 crashing, is unusable

[D. S. McNeil]

 For future reference, run `md5sum filename` in a terminal to check the
 MD5 sum of a file.

I'm pretty sure it's md5, not md5sum, on OS X (at least in 10.6).


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Sage 5.0 crashing, is unusable

[D. S. McNeil]

I can reproduce this on my 10.6.8 macbook:

sage: int(2)
2
sage: int(2.75)

Program received signal EXC_BAD_INSTRUCTION, Illegal instruction/operand.
0x000101723ed9 in case1 ()
(gdb) bt
#0  0x000101723ed9 in case1 ()
#1  0x000103e8bba4 in parsed_string_to_mpfr ()
#2  0x000103e8c7fb in mpfr_strtofr ()
Previous frame inner to this frame (gdb could not unwind past this frame)

Since this is the computer I actually use, I keep things in /usr/local
and /opt/local (including copies of gmp), so I tried moving those in
case it was picking up something it shouldn't, but no luck.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Sage 5.0 crashing, is unusable

[D. S. McNeil]

 Doug, are you using the binary or did you compile your Sage?

Deliberately the binary; I've never had problems with a Sage I've
successfully compiled myself.   [Haven't compiled 5.0 myself yet on
the Mac, though I did at work today on ubuntu 12.04 and it went fine.
Will probably try overnight.]

This feels sort of like a 32-bit/64-bit issue, and I'm getting some
unexpected results in some directions.  Trying to remind myself about
some of the Mac quirks.

FWIW even sage: 2.3 trips the error for me, so I think pretty much
every mpfr library call will cause troubles.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Sage 5.0 crashing, is unusable

[D. S. McNeil]

I was able to find a case which crashed but doesn't have to crash
Sage.  Doesn't tell us much new, but here it is anyway:


sage: 12345678912345678912345678912345678901234567890.
---
RuntimeError  Traceback (most recent call last)

/Users/mcneil/Desktop/ipython console in module()

/Applications/Sage-5.0-OSX-64bit-10.6.app/Contents/Resources/sage/local/lib/python2.7/site-packages/IPython/Prompts.pyc
in __call__(self, arg)
550
551 # and now call a possibly user-defined print mechanism
-- 552 manipulated_val = self.display(arg)
553
554 # user display hooks can change the variable to be stored in

/Applications/Sage-5.0-OSX-64bit-10.6.app/Contents/Resources/sage/local/lib/python2.7/site-packages/IPython/Prompts.pyc
in _display(self, arg)
576 return IPython.generics.result_display(arg)
577 except TryNext:
-- 578 return self.shell.hooks.result_display(arg)
579
580 # Assign the default display method:

/Applications/Sage-5.0-OSX-64bit-10.6.app/Contents/Resources/sage/local/lib/python2.7/site-packages/IPython/hooks.pyc
in __call__(self, *args, **kw)
139 #print prio,prio,cmd,cmd #dbg
140 try:
-- 141 ret = cmd(*args, **kw)
142 return ret
143 except ipapi.TryNext, exc:

/Applications/Sage-5.0-OSX-64bit-10.6.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/misc/displayhook.pyc
in result_display(ip_self, obj)
150 # IPython's default result_display() uses the
IPython.genutils.Term.cout stream.
151 # See also local/lib/python2.6/site-packages/IPython/hooks.py.
-- 152 print_obj(IPython.genutils.Term.cout, obj)
153
154 def displayhook(obj):

/Applications/Sage-5.0-OSX-64bit-10.6.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/misc/displayhook.pyc
in print_obj(out_stream, obj)
142 if _check_tall_list_and_print(out_stream, obj):
143 return
-- 144 print out_stream, `obj`
145
146 def result_display(ip_self, obj):

/Applications/Sage-5.0-OSX-64bit-10.6.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/rings/real_mpfr.so
in sage.rings.real_mpfr.RealNumber.__repr__
(sage/rings/real_mpfr.c:9440)()

/Applications/Sage-5.0-OSX-64bit-10.6.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/rings/real_mpfr.so
in sage.rings.real_mpfr.RealNumber.str
(sage/rings/real_mpfr.c:11259)()

RuntimeError: Illegal instruction

or under gdb:


sage: 12345678912345678912345678912345678901234567890.

Program received signal EXC_BAD_INSTRUCTION, Illegal instruction/operand.
0x000101723ec7 in case3 ()
(gdb) bt
#0  0x000101723ec7 in case3 ()
#1  0x000103e6554f in mpfr_get_str ()
#2  0x000103dfac65 in __pyx_pf_4sage_5rings_9real_mpfr_10RealNumber_13str ()


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-notebook] Re: [sage-support] Re: Exporting data from sage notebook

[D. S. McNeil]

 I think you should close that.  The following works fine in the
 notebook.   It makes no sense to close the file -- instead, you have
 to delete the csv writer object, which flushes it to the file.

This isn't guaranteed to flush it to the file, though, because del
only deletes the name.  That it results in a flush is entirely a
cpython implementation detail because there are no live references to
it.  Looking at the test.csv file after `del output`:

mcneil@ubuntu:~/coding$ python csvflush.py
['1,2\r\n', '3,4\r\n', '5,6\r\n']
mcneil@ubuntu:~/coding$ pypy csvflush.py
[]

Or more specifically in cpython:

import csv
fp = open(test.csv, w)
output=csv.writer(fp)
data=[[1,2],[3,4],[5,6]]
for r in data:
   output.writerow(r)
print(open('test.csv').readlines())
del output
print(open('test.csv').readlines())
del fp
print(open('test.csv').readlines())

produces

[]
[]
['1,2\r\n', '3,4\r\n', '5,6\r\n']

where nothing happens until the last reference to the *file* is lost,
which in cpython gets it closed.

I'd either call .flush() and/or .close() explicitly or wrap it in a
with block so that I don't have to remember.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Bug in edge_cut of undirected weighted graphs

[D. S. McNeil]

I think I have a smaller example:

sage: G = Graph([(0, 3, 1), (0, 4, 1), (1, 2, 1), (2, 3, 1), (2, 4, 1)])
sage: G.edge_cut(0,1,value_only=False,use_edge_labels=True)
[1, [(0, 3, 1), (1, 2, 1), (2, 3, 1)]]
sage: G.edge_cut(0,1,value_only=False,use_edge_labels=True,method='LP')
(1.0, [(1, 2)])

This should make debugging easier.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] About the solve function

[D. S. McNeil]

You're probably getting a list of dictionaries:

sage: var(x y)
(x, y)
sage: sols = solve([x^2-4==0, y^2-9==0], x, y, solution_dict=True)
sage: sols
[{y: -3, x: -2}, {y: -3, x: 2}, {y: 3, x: -2}, {y: 3, x: 2}]

This is just like any other list -- the fact the elements happen to be
dictionaries doesn't change anything -- and can be looped over and
indexed.  For example:

sage: sols[0]
{y: -3, x: -2}
sage: sols[0][x]
-2
sage: sols[0][y]
-3
sage: sols[1][x]
2
sage: sols[1][y]
-3
sage: sum(sol[x]+sol[y] for sol in sols)
0

etc.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Legend label problem

[D. S. McNeil]

This seems to work on test.sagenb.org (5.0.beta1) and the 5.0.beta4 I
have around, so something (whether Sage-side or matplotlib-side) must
have changed for the better from 4.8.  Is it time-sensitive enough to
track it down and make a backpatch?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Question about plot

[D. S. McNeil]

 I am a two-tuple vector, vet=[(1,2),(3,4),(5,6),..], i want plot
 this data with a line aproximation, (interpolation this points), exist
 any parameter in list_plot function for this.

Maybe

line(vet)

does what you want?  You can look at
http://www.sagemath.org/doc/reference/sage/plot/plot.html for other
possibilities.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] gcd's of numbers mod N

[D. S. McNeil]

 sage seems to think that the gcd of 6 and (-2 mod 6) is -2 mod 6, which it 
 converts to 4.  A mathematician would say that the gcd is 2.
 Is this a bug, or does sage have a higher purpose here?

Sage is actually reasoning slightly differently, I think.  First it
decides whether there's a canonical coercion to a common parent.  In
this case, it concludes that the common parent should be the ring of
integers modulo 6:

sage: parent(Mod(4,6))
Ring of integers modulo 6
sage: parent(6)
Integer Ring
sage: z = cm.canonical_coercion(Mod(4,6), 6)
sage: z
(4, 0)
sage: parent(z[0]), parent(z[1])
(Ring of integers modulo 6, Ring of integers modulo 6)

There's no gcd method defined in this ring, so it falls back to
attempting to coerce 4 mod 6 and 0 mod 6 to ZZ, which succeeds, and we
get the integer version under which we have

sage: gcd(4,0)
4

It's not clear to me what the best way to handle this case is.  Paging
Simon King.. :^)


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: speed question, numpy vs CDF

[D. S. McNeil]

Anyone using numpy from Sage should beware of the following:

sage: import numpy
sage: m = numpy.matrix([[1,2],[3,4]])
sage: m[:,0]
matrix([[1, 3]])
sage: m[:,int(0)]
matrix([[1],
[3]])

That is, if you use a Sage integer to index a numpy matrix, you don't
get the expected shape back.  I know exactly why this happens --
http://trac.sagemath.org/sage_trac/ticket/10928 -- but ideally it
should be patched upstream.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] sum() - TypeError: cannot evaluate symbolic expression numerically

[D. S. McNeil]

On Sun, Feb 12, 2012 at 2:13 PM, mrkvon michal.sala...@gmail.com wrote:
 Hello,

 what is wrong with this, please? I tried to sum some matrix elements
 and got error below. I tried to sum by other means (comment # and ##
 in example below) without problem. I even defined my own summation
 function (analogic to ##), which at the end made similar error.
 I don't think I try to evaluate numerically, so why this error?

Actually, you _are_ trying to evaluate numerically, simply not what
you thought you were. :^)

The problem with

sum(lore4[k,k],k,0,3)

isn't the sum; it's the lore4[k,k]:


sage: lore4[k,k]
[...]
TypeError: cannot evaluate symbolic expression numerically
---

It's trying to convert the k to an integer index to access the matrix
elements.  By the way, whenever Sage (or Python) gives me a TypeError
I break everything down into its pieces like this and try to find
which one is causing the problems.  It doesn't always work, but I'd
say 75% of the time it's an input formatting error on my part that I
can find.

Anyway, I'm not sure if Sage matrices can be indexed by symbolic
variables like this.  [Whenever I say you can't do something in Sage
someone proves me wrong, but I'm not familiar with the trick if there
is one.]

In this case, rather than using the symbolic sum 4-term calling style
(expression, var, v0, v1), I'd write it with a generator expression
instead:

---
sage: sum(lore4[k,k] for k in [0..3])
2*sinh(x) + 2
---

or something like that.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] finding a p-cycle of a real function

[D. S. McNeil]

Warning: I haven't thought this through, but the ideas might be useful.

It looks like you'd call something a 3-cycle if for some x_0 we had
f(f(f(x_0))) == x_0, right?  Then we should be able to do this
numerically, with some caveats:

# can't remember the slick way, so brute force
def iter_apply(f0, n):
fs = [f0()]
for i in xrange(n-1):
last = fs[-1]
fs.append(f0(last))
return fs

def find_pcycles(f0, n):
fs = iter_apply(f0, n)
req = fs[-1] == x # defining equation of the cycle
roots = req.roots(ring=RR)
for root, mult in roots:
yield [fi(x=root) for fi in fs]

which gives

sage: f(x) = x^2-1.76
sage:
sage: p1 = list(find_pcycles(f, 1))
sage: p1
[[-0.917744687875782], [1.91774468787578]]
sage: f(p1[0][0])
-0.917744687875783
sage: f(p1[1][0])
1.91774468787578
sage: list(find_pcycles(f, 2))
[[0.504987562112089, -1.50498756211209], [-0.917744687875783,
-0.917744687875782], [-1.50498756211209, 0.504987562112089],
[1.91774468787578, 1.91774468787578]]
sage: list(find_pcycles(f, 3))
[[1.33560128916886, 0.0238308036295167, -1.75943209279837],
[1.27545967679486, -0.133202612870350, -1.74225706392451],
[-0.917744687875783, -0.917744687875782, -0.917744687875783],
[-1.74225706392451, 1.27545967679486, -0.133202612870348],
[-1.75943209279837, 1.33560128916886, 0.0238308036295143],
[-0.133202612870349, -1.74225706392451, 1.27545967679486],
[0.0238308036295150, -1.75943209279837, 1.33560128916886],
[1.91774468787578, 1.91774468787578, 1.91774468787578]]

Note that this says [-0.917744687875783, -0.917744687875782,
-0.917744687875783] is a 3-cycle, which it kind of is..

Anyway, maybe playing off of something like that would help?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Numerical approximation of symbolic coefficients

[D. S. McNeil]

 First attempt: loop through each term and try to n() the coefficient.
 Madness.

Based on a suggestion Mike Hansen once gave me --
http://ask.sagemath.org/question/411/substituting-expressions-for-numbers
-- I tend to use subclasses of Converter when I need to do something
like this, so as not to get lost in the madness. :^)  Something like:


from sage.symbolic.expression_conversions import Converter

class Evaluator(Converter):
def arithmetic(self, ex, operator):
return reduce(operator, map(self, ex.operands()))
def pyobject(self, ex, obj):
return ex.n()

sage: E = Evaluator()
sage: var(u x)
(u, x)
sage: q = ((314069483520)*sqrt(3/(sin(u+2)))*u - 80295755776*x +
4831838208)/(1953125*x^63)
sage: q
33554432/1953125*(9360*sqrt(3)*u*sqrt(1/sin(u + 2)) - 2393*x + 144)/x^63
sage: E(q)
17.17986918398*(16211.99555884469*u*(1/sin(u + 2.0))^0.5 -
2393.0*x + 144.0)/x^63.0



Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] issue with combining sub lists into a list

[D. S. McNeil]

 [[1],[2,8,4,6],[1,2,7],[9,3,4,6],...]
 [[1],[14,17,18,19],[1,4,11],[9,14,16,19],...]

 would like to figure out how to get both results like this:

 [1,2,8,4,6,1,2,7,9,3,4,6,...]
 [1,14,17,18,19,1,4,11,9,14,16,19...]

Two ways come to mind:

sage: a = [[1],[2,3,4],[5,6,7]]
sage: a
[[1], [2, 3, 4], [5, 6, 7]]
sage: flatten(a)
[1, 2, 3, 4, 5, 6, 7]
sage: sum(a,[])
[1, 2, 3, 4, 5, 6, 7]

The two approaches differ in how deep they go: flatten recurses
downwards whereas the sum only goes down one level.  Sometimes you
want one or the other:

sage: b = [[1],[[2,3],[4,5]]]
sage: b
[[1], [[2, 3], [4, 5]]]
sage: flatten(b)
[1, 2, 3, 4, 5]
sage: sum(b,[])
[1, [2, 3], [4, 5]]


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: help with filtering between two lists.

[D. S. McNeil]

 Most of the time, filtering jobs like this can be accomplished by one-
 liners using the right list combination and iteration tools:

 [c for c in enumerate(zip(l1,l2)) if c[1][0] ==c[1][1]]

You can even avoid the [1][0] stuff by using some nifty unpacking:

[(c,ai) for c, (ai,bi) in enumerate(zip(a,b)) if ai == bi]


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Plotting weighted graph

[D. S. McNeil]

 Hi all, I'm confused why this code results in a graph plotted with the
 (single) edges labeled with the weight between each pair of vertices:

 M = Matrix([[0,1,-1,5],[1,0,-1/2,-1],[-1,-1/2,0,2],[5,-1,2,0]])
 G = Graph(M,sparse=True)
 G.plot(edge_labels=True)

 Whereas this code does not and instead has multi-edges labeled None:

 M = Matrix([[0, 1, 0, 0], [1, 0, 1, 2], [0, 1, 0, 0], [0, 2, 0, 0]])
 G = Graph(M,sparse=True)
 G.plot(edge_labels=True)

 What am I missing?

I think it's something like this (take with grain of salt, I'm bad at
multitasking):

Briefly, it's not interpreting the second matrix as a weighted
adjacency matrix.  If you specifically tell it that it is, say by

M = Matrix([[0, 1, 0, 0], [1, 0, 1, 2], [0, 1, 0, 0], [0, 2, 0, 0]])
G = Graph(M,weighted=True)

it should behave.

Less briefly: the code in graph.py uses the following logic to decide
whether these two are weighted or not:

if format == 'adjacency_matrix':
entries = uniq(data.list())
for e in entries:
try:
e = int(e)
assert e = 0
except:
if weighted is False:
raise ValueError(Non-weighted graph's+
 adjacency matrix must have only nonnegative+
 integer entries)
weighted = True
if multiedges is None: multiedges = False
break

IOW, to decide if it should set weighted to True, it loops over each
entry and (1) tries to convert it to an integer, and (2) asserts that
it's non-negative.  If either of those fail, then it decides that
weighted should be True (in this case weighted starts as None).

So the first graph works because the nonnegative check fails (or the
integral check fails, I'm too lazy to confirm what direction the
matrix entries are iterated in..)  The second graph's matrix is
integral and nonnegative, so it sneaks through, and things get silly
after that.

Would you like to open a bug report on trac?  This should be easy to fix.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] sagenb down

[D. S. McNeil]

test.sagenb.org seems to be up, so if she has a saved local copy of a
worksheet, she could use that..


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] matrix help

[D. S. McNeil]

 def U(N,M):
    U=matrix(ZZ,N*M)
    for i in range(N*M):
        for j in range(N*M):
            U[i,j]=1
    return U
 def Q(N,M):
    Q=matrix(ZZ,N*M)
    for i in range(N*M):
        for j in range(N*M):
            Q[i,j]=U(N,M)
    return Q

U(N,M) is a function which returns a matrix.  The line Q[i,j]=U(N,M)
attempts to set the (i,j)-th entry of the Q matrix, defined over ZZ,
to an NxM matrix, which isn't going to work.  (Hence the cannot
coerce a matrix to an integer error message -- it's trying to fit
U(N,M) into an entry of Q, and it can't figure out how to do it.)

BTW, you might also be interested in the command ones_matrix:

sage: ones_matrix(ZZ,3,4)
[1 1 1 1]
[1 1 1 1]
[1 1 1 1]


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Why this difference with mpmath numerical differentiation?

[D. S. McNeil]

 Normally I can differentiate analytically to get the slope and
 inflection points. While trying to have a quick look into these points
 with numerical differentiation I noticed that mpmath is giving me
 values quite distant from the ones obtained with analytical formulas.

I think there's a missing term in your analytical derivatives.  Your
three functions are

# fundamental mode
def phi(x):
   return mpmath.cos(k*x) - mpmath.cosh(k*x) + alpha* (
mpmath.sin(k*x) - mpmath.sinh(k*x) )

# slope of fundamental
def phi_dx(x):
   return -mpmath.sin(k*x) - mpmath.sinh(k*x) + alpha* (
mpmath.cos(k*x) - mpmath.cosh(k*x) )

# curvature of fundamental
def phi_dx2(x):
   return -mpmath.cos(k*x) - mpmath.cosh(k*x) - alpha* (
mpmath.sin(k*x) + mpmath.sinh(k*x) )

and it looks to me like you're missing a factor of k in phi_dx and k^2
in phi_dx2.  Inserting them seems to get rid of the discrepancy.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Hole in algorhitm

[D. S. McNeil]

:-/ That's definitely a bug.  The precision used in
integral_points_with_bounded_mw_coeffs (100 bits) is too small to find
that solution.  Even bumping the precision to 120 bits suffices,
although that's probably the wrong approach:

[(436 : 559 : 1), (450 : 2925 : 1), (666 : 14589 : 1), (900 : 25425 :
1), (1150 : 37925 : 1), (1800 : 75825 : 1), (2619 : 133722 : 1),
(26154 : 4229667 : 1), (27675 : 4603950 : 1), (748476100 :
20477027135825 : 1)]

This is now http://trac.sagemath.org/sage_trac/ticket/12095


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] 'unable to simplify to float approximation' error in Logistic equation

[D. S. McNeil]

The problem is that the solution that desolve returns:

 sage: myode = tau*diff(p,t) == p*(1-p/k)
 sage: sol(t) = desolve(de=myode, ivar=t, dvar=p)

 sage returns:
  -tau*log(-k + p(t)) + tau*log(p(t)) == c + t

isn't in an easy enough form for plot to display.  You're setting
sol(t) not to an expression but to an equation, and a non-trivial one.
 (You might be able to make implicit_plot work with some tweaking but
I'm lazy. :^)

We can fix that, though.  Try something like



tau = var('tau') # proportionality constant
k = var('k') # max population size
t = var('t') # independent variable
p = function('p', t) # dependent variable

myode = tau*diff(p,t) == p*(1-p/k)
sol = desolve(de=myode, ivar=t, dvar=p)
sol = sol.simplify_full() # simplify_log would suffice
sol_p(t) = solve(sol,p,solution_dict=True)[0][p] # make a function
from the solution
plot(sol_p.subs(tau=20, k=10, c=1), 0, 100)




Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Error on variables used by find_fit

[D. S. McNeil]

You're assuming that the fit function is going to return the variables
in the same order in both cases, so that (e.g.) a and k1 will be the
first in both cases.  That's not true, though.  I found:

[a == 3.4102956225507519, b == 0.010344876276231638, c ==
-0.00094076216744204172, d == (1.8097013599216266e-05), f ==
(-1.7547848824431721e-07), g == (1.0168442135198695e-09), h ==
(-3.7548063250237949e-12), j == (8.9493796546766707e-15), k ==
(-1.3366803967120946e-17), l == (1.1398034585075684e-20), n ==
(-4.2382572237693009e-24), o == 3.0511588067264963]

[b == 0.010358179948437265, c == -0.00094126295808369734, d ==
(1.8106258279315326e-05), f == (-1.7557568485259962e-07), g ==
(1.0174721357045694e-09), h == (-3.7573869398464435e-12), j ==
(8.9561476716220824e-15), k == (-1.331928232302e-17), k1 ==
3.4101714386786477, l == (1.1408041039846924e-20), n ==
(-4.242187409468519e-24), o == 3.0512548011942671]

so the indexing is all off.  Not surprising the second one explodes.
Return order looks to be alphabetical, but I'm too lazy to find out if
that's guaranteed.  In any case, we can use a solution_dict instead to
make life much easier.

So what I'd recommend, given your V, is something like this:

sage: N_order = 12
sage:
sage: # make a list of symbolic variables
sage: coeffs = list(var(c_%d % i) for i in range(N_order))
sage: coeffs
[c_0, c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8, c_9, c_10, c_11]
sage:
sage: # build the model [here just a polynomial]
sage: model(x) = sum(c_i*x**i for i, c_i in enumerate(coeffs))
sage: model
x |-- c_11*x^11 + c_10*x^10 + c_9*x^9 + c_8*x^8 + c_7*x^7 + c_6*x^6 +
c_5*x^5 + c_4*x^4 + c_3*x^3 + c_2*x^2 + c_1*x + c_0
sage:
sage: # actually fit it, returning a dictionary
sage: mfit = find_fit(V, model, solution_dict=True)
sage: mfit
{c_7: -3.7564661583891436e-12, c_0: 3.0514755095206816, c_6:
1.0172525553243703e-09, c_11: -4.2407560895324755e-24, c_10:
1.1404399242138324e-20, c_3: -0.00094112724266502126, c_9:
-1.3373790954916053e-17, c_2: 0.010355836374653954, c_4:
1.8103350226003258e-05, c_1: 3.4101603849319178, c_8:
8.9537056388588457e-15, c_5: -1.7554291229990257e-07}
sage:
sage: # substitute these values into the model, getting a function
sage: mfunc = model.subs(mfit)
sage: mfunc
x |-- -(4.2407560895324755e-24)*x^11 + (1.1404399242138324e-20)*x^10
- (1.3373790954916053e-17)*x^9 + (8.9537056388588457e-15)*x^8 -
(3.7564661583891436e-12)*x^7 + (1.0172525553243703e-09)*x^6 -
(1.7554291229990257e-07)*x^5 + (1.8103350226003258e-05)*x^4 -
0.00094112724266502126*x^3 + 0.010355836374653954*x^2 +
3.4101603849319178*x + 3.0514755095206816
sage:

This way I don't have to worry about keeping all those letters
straight.  Note though that unless you're sure it really is a
polynomial, in general fitting data with a high-order polynomial is a
very bad idea (pretty unstable).

Does that help?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: ValueError when using np.random.seed()

[D. S. McNeil]

thread hijack but not really

 Most likely, the problem here is that numpy doesn't recognize Sage Integer
 objects as being convertable to python integers.  I think this is a failing
 of numpy (it should check the __index__ method to see if it can convert the
 Integer to a python integer).

Speaking of numpy issues, does someone with more influence than I have
(i.e. anyone with nonvanishing influence) want to reraise the
http://trac.sagemath.org/sage_trac/ticket/10928 issue with the numpy
guys?

/thread hijack


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Letter to Number mapping

[D. S. McNeil]

Do you know about Python dictionaries?

http://docs.python.org/tutorial/datastructures.html#dictionaries

They're a kind of general map along the lines that you want, and you
can build them in many ways.  For example:

import string
m = dict((c, ord(c.lower())-ord('a')+17) for c in string.ascii_letters)
print m['A'], m['a']

or

indices = range(5, 26+5)
m2 = dict(zip(string.ascii_lowercase, indices))
print m2['a'], m2['z']

or you can build one manually:

m3 = {'a': 5, 'b': 238}
print m3['a'], m3['b']


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Bessel Contour Plot problem

[D. S. McNeil]

Hi!

I'm not sure what you mean by

 Also if I have the terms, and return on separate lines in the
 definition then I get this error.

You _have_ to have the different statements on separate lines, like I
did.  (Unless you use a semicolon to separate the statements, I
suppose.)  Did you combine them?  It looks like you did but it's hard
to be sure because sometimes carriage returns are lost during the
cut-and-paste process.

Paying attention to the syntax is even more important Python than in
many other languages because Python is whitespace-sensitive-- if you
don't use the right indentation, the code might not work.  I think it
might be useful to read through a Python tutorial if you're more
familiar with the syntax from other languages.  It won't take too long
but may clear up some of the confusion.

As for the NameError:

 Well I did use your method but I still get an error.

Actually, it doesn't look like you did use my code.  You modified the
last line to read return term1 + term2 if t != 0 else infinity,
thereby introducing a variable t which you hadn't defined, which leads
to the NameError explaining -- reasonably enough -- that you're
referring to a name 't' which you didn't define.

Does that make sense?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Bessel Contour Plot problem

[D. S. McNeil]

 def f(y,t): return (1+(bessel_J(0, gro)/bessel_J(2, gro)))*(r/
 kro)*(bessel_J(1, r)/bessel_J(1,gro))*z.cos()-(bessel_J(0,r)/
 bessel_J(2, gro))*(r**2/kro**2) if t != 0 else infinity

A few things:

(1) Your arguments to this function are y and t, but inside you use r
and z.  r and z are still the symbolic values.  That's why the error
message says Cannot evaluate symbolic expression to a numeric value,
and that's usually the problem: some variable isn't being assigned a
number.

(2) Since z is going to be a float, not a Sage real, z.cos() isn't
going to work as .cos() isn't a float method.  (Sorry about this; in a
sense it's an implementation detail there's no way to guess.  Python
floats are low-level dumb numbers which only have a few methods
living inside them.)

(3) The if condition else infinity trick was a workaround for a
singularity that wasn't being handled too well.  If there isn't one,
there's no need to worry about it.

So I'd simply replace your f function with something like

def f(r,z):
term1 = (1+(bessel_J(0, gro)/bessel_J(2,
gro)))*(r/kro)*(bessel_J(1, r)/bessel_J(1,gro))*cos(z)
term2 = -(bessel_J(0,r)/bessel_J(2, gro))*(r**2/kro**2)
return term1 + term2

[possible typos, I didn't check too carefully, but you get the idea]


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Bessel Contour Plot problem

[D. S. McNeil]

 Sage: y,t = var('y,t')
 Sage: contour_plot(lambda y,t: (sqrt(t^2+y^2)/(2*pi*y))*(bessel_J(0,
 t ).arccos()), (t, 0, 3), (y, 0, 4), fill = false, axes_labels=['$gro
 $','$kro$'], contours = [0.1, 0.2, 0.3, 0.5, 0.7, 1.0, 2.0, 4.0],
 fill=False, legend_label='qwall');

 This is the error message that I get when processing the code:

 The error I get is: File ipython console, line 1 SyntaxError:
 keyword argument repeated (ipython console, line 1).

 Does it have something to do with assigning the variables?

It has to do with the fact you're repeating a keyword argument. :^)
Note that you write fill=false and then later fill=False.

Also, be careful: when you use a Python function (lambda or def), you
can't use the (variable, lower, upper) syntax, you have to give them
in order.  Proof:

sage: aa = []
sage: bb = []
sage: def f(a,b):
: aa.append(a)
: bb.append(b)
: return a+b
:
sage: a,b = var(a b)
sage: contour_plot(f,(b,2,3),(a,4,5))

sage:
sage: min(aa),max(aa)
(2.0, 3.0)
sage: min(bb),max(bb)
(4.0, 5.0)

And you see that f gets two variables ranging from 2 to 3 and 4 to 5,
but it doesn't put the a values in a.  (Aside for devs: we could fix
this using introspection.)

Finally, I found that I had to manually handle the zero case to avoid
getting a Symbolic division by zero error:


y,t = var('y,t')

def f(y,t):
return (sqrt(t^2+y^2)/(2*pi*y))*(bessel_J(0,t).arccos()) if y != 0
else infinity

p = contour_plot(f, (y, 0, 6), (t, 0, 3), fill = false,
axes_labels=['$gro$','$kro$'], contours = [0.1, 0.2, 0.3, 0.5, 0.7,
1.0, 2.0, 4.0], legend_label='qwall');
show(p)


Hope some of that helps!


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Question about legend place in plot

[D. S. McNeil]

 I want that legend is in left side, ... how?

Try something like:

p = plot(sin,legend_label=sin)
p.set_legend_options(loc='upper left')
p.show()

After making some plot p, you can type help(p.set_legend_options) to
see more information about the various things you can configure.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Question about Patterson Algorithm Implementation

[D. S. McNeil]

 This is definitely not a bug.   The definition of the _add_ method
 absolutely demands that both inputs have exactly the same parent.  In
 the above instance, the left hand input (=1) has parent ZZ, and the
 right hand input (=SR(2)) has parent the symbolic ring.

Yeah, I know that-- it's the violation of that assumption which
ultimately crashed the OP's code, after all.

I guess I've inherited the bias from Python that users shouldn't be
able to segfault the interpreter from pure Python code.
Anything Cythonic probably falls into the Sage equivalent of the
ctypes exception class, and I guess you can get the same crash with
any non-typechecking cpdef'd object, but it still feels wrong.

Meh.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Question about Patterson Algorithm Implementation

[D. S. McNeil]

I don't think you should need to call _add_, but this looks like a bug to me:

--
| Sage Version 4.7.1, Release Date: 2011-08-11   |
| Type notebook() for the GUI, and license() for information.|
--
sage: 1+SR(2)
3
sage: 1.__add__(SR(2))
3
sage: 1._add_(SR(2))


Unhandled SIGSEGV: A segmentation fault occurred in Sage.
This probably occurred because a *compiled* component of Sage has a bug
in it and is not properly wrapped with sig_on(), sig_off(). You might
want to run Sage under gdb with 'sage -gdb' to debug this.
Sage will now terminate.



Similarly for _sub_ and _mul_.  _div_ gives a slightly different result:


sage: 1._div_(SR(2))
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (3258, 0))

ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (3244, 0))

---
TypeError Traceback (most recent call last)

/Users/mcneil/ipython console in module()

/Applications/sage/local/lib/python2.6/site-packages/sage/rings/integer.so
in sage.rings.integer.Integer._div_ (sage/rings/integer.c:11866)()

/Applications/sage/local/lib/python2.6/site-packages/sage/rings/integer.so
in sage.rings.integer.Integer._div_ (sage/rings/integer.c:11824)()

/Applications/sage/local/lib/python2.6/site-packages/sage/rings/integer_ring.so
in sage.rings.integer_ring.IntegerRing_class._div
(sage/rings/integer_ring.c:5158)()

TypeError: Argument 'right' has incorrect type (expected
sage.rings.integer.Integer, got sage.symbolic.expression.Expression)


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: error: C preprocessor /lib/cpp fails sanity check (installing sage-4.7.1 from source)

[D. S. McNeil]

 The shortest way to test it my above interpretation of your logs is indeed
 correct is by seeing if excecuting:
 echo #include limits.h | cpp

 gives any errors. If so my interpretation is correct.

And if it is, which seems likely, there's probably a kernel headers
package missing.  Googling suggests it's probably called
kernel-headers in Pardus, but maybe kernel-module-headers or
kernel-source or something.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Chinese Remainder Theorem

[D. S. McNeil]

On Fri, Sep 23, 2011 at 12:39 AM, Santanu Sarkar
sarkar.santanu@gmail.com wrote:
 I want to find integer such that
 x= 1 mod 3
 x=2  mod 5
 x=3  mod 7
 like this system of congruences using Chinese Remainder Theorem.
 In Sage, crt() function takes only 4 argument.

sage: help(CRT)

crt(a, b, m=None, n=None)
Returns a solution to a Chinese Remainder Theorem problem.

INPUT:

- ``a``, ``b`` - two residues (elements of some ring for which
  extended gcd is available), or two lists, one of residues and
  one of moduli.
[...]

If ``a`` and ``b`` are lists, returns a simultaneous solution to
the congruences `x\equiv a_i\pmod{b_i}`, if one exists.

.. SEEALSO::

- :func:`CRT_list`


sage: CRT([1,2,3],[3,5,7])
52
sage: x = CRT([1,2,3],[3,5,7])
sage: x % 3, x % 5, x % 7
(1, 2, 3)


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] How to write Sage code to cython code

[D. S. McNeil]

 It always returns 101, not a random prime of 100 bit integer.

That's because in Python/Cython, the carat ^ isn't exponentiation,
it's bitwise xor.  The most general solution is to use **:

Python 2.7.2 (v2.7.2:8527427914a2, Jun 11 2011, 15:22:34)
[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] on darwin
Type help, copyright, credits or license for more information.
 2^(100-1)
97
 2^(100)
102
 2**(100-1)
633825300114114700748351602688L
 2**(100)
1267650600228229401496703205376L


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Obtaining Size of a List and of elements of a List

[D. S. McNeil]

In your code, ComSet is a Python list (not a set) as are many of its
components, and you use len(x) to get the size:

sage: ComSet, type(ComSet), len(ComSet)
([[[0, 1], [0, 2], [1, 2]], [[0, 1, 2]], [[0, 1], [0, 2], [1, 2]]],
type 'list', 3)
sage: ComSet[0], type(ComSet[0]), len(ComSet[0])
([[0, 1], [0, 2], [1, 2]], type 'list', 3)
sage: ComSet[0][0], type(ComSet[0][0]), len(ComSet[0][0])
([0, 1], type 'list', 2)
sage: ComSet[0][0][0], type(ComSet[0][0][0])
(0, type 'int')

cardinality is a method not of Python lists, but of the Combinations
object.  For example:

sage: C
Combinations of [0, 1, 2] of length 2
sage: C.cardinality()
3
sage: list(C)
[[0, 1], [0, 2], [1, 2]]
sage: C.list()
[[0, 1], [0, 2], [1, 2]]
sage: len(C.list())
3

The reason tab-completion doesn't reveal len is because len is a
function, not a method on the object, and the dot-tab procedure
returns the object's contents.  (Admittedly, if you type
ComSet.__[tab], you can see the special methods, including
ComSet.__len__ which is used behind the scenes, but you would never
write ComSet.__len__() in real code.)

Does that help?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Sums of Modular Symbols

[D. S. McNeil]

 Thanks for this.  There still seems to be a manual step in going from,
 say,
 s1 = 2*(1,8) - (1,9)
 to
  s1 = 2*b[1] - b[2]

I may be misunderstanding you.  Are you saying you want to enter the line

s1 = 2*(1,8)-(1,9)

verbatim and have it work?  That I don't think I can do (unless you're
willing to type them not at the prompt but in a string or file which
is then parsed).  One option is to use a dictionary as an intermediary
step:

age: M = ModularSymbols(11,2);M
Modular Symbols space of dimension 3 for Gamma_0(11) of weight 2 with
sign 0 over Rational Field
sage: m = dict((sage_eval(str(b)), b) for b in M.basis())
sage: m
{(1, 0): (1,0), (1, 8): (1,8), (1, 9): (1,9)}

This dictionary basically maps a tuple to the corresponding
ModularSymbolsElement in the basis.  [Yes, to be safer I should really
use ast.literal_eval.]  Then you could write:

sage: s = [2*m[1,8] - m[1,9],
:  -m[(1,0)] + m[(1,9)],
:  -m[(1,0)] + m[(1,8)]]
sage: s
[2*(1,8) - (1,9), -(1,0) + (1,9), -(1,0) + (1,8)]

Note that the parentheses () are optional.  Finally:

sage: V = span([sx.list() for sx in s],QQ)
sage: V
Vector space of degree 3 and dimension 3 over Rational Field
Basis matrix:
[1 0 0]
[0 1 0]
[0 0 1]


Would something like that work or do you have a particular format that
you can't escape?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Sums of Modular Symbols

[D. S. McNeil]

.. I suppose you could even add a convenience function

def b(*x):
return m[x]

after which

sage: s = [2*b(1,8) - b(1,9),
:  -b(1,0) + b(1,9),
:  -b(1,0) + b(1,8)]
sage: s
[2*(1,8) - (1,9), -(1,0) + (1,9), -(1,0) + (1,8)]

would work.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Sums of Modular Symbols

[D. S. McNeil]

 We would like to know if certain sums of modular symbols span the
 space.

Is this the sort of thing you had in mind?

sage: M=ModularSymbols(11,2);M
Modular Symbols space of dimension 3 for Gamma_0(11) of weight 2 with
sign 0 over Rational Field
sage: b = M.basis()
sage:
sage: s1 = 2*b[1] - b[2]
sage: s2 = -b[0] + b[2]
sage: s3 = -b[0] + b[1]
sage: s1, s2, s3
(2*(1,8) - (1,9), -(1,0) + (1,9), -(1,0) + (1,8))
sage: s1.list()
[0, 2, -1]
sage:
sage: V = span([s.list() for s in s1, s2, s3],QQ); V
Vector space of degree 3 and dimension 3 over Rational Field
Basis matrix:
[1 0 0]
[0 1 0]
[0 0 1]
sage:
sage: s3 = s2
sage: V = span([s.list() for s in s1, s2, s3],QQ); V
Vector space of degree 3 and dimension 2 over Rational Field
Basis matrix:
[   10   -1]
[   01 -1/2]


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Nearest Integer

[D. S. McNeil]

 How to find the nearest integer (+ve or -ve) of a rational number (P/Q)
 where P,Q are very large integers?

You could use the .round method of rationals.

sage: q = 17+1/2+1/11**1000
sage: RR(q.numerator()), RR(q.denominator())
(2.48685403212345e10413928, 1.42105944692768e10413927)
sage: q.round()
18

sage: q = 13**(10**7)+1/2+1/11**(10**7)
sage: RR(q.numerator()), RR(q.denominator())
(4.73893349282680e21553360, 1.42105944692768e10413927)
sage: z = q.round()
sage: factor(z-1)
13^1000


sage: help(q.round)
round(...)
File: sage/rings/rational.pyx (starting at line 2858)

Returns the nearest integer to self, rounding away from 0 by
default, for consistency with the builtin Python round.

INPUT:


-  ``self`` - a rational number

-  ``mode`` - a rounding mode for half integers:

   - 'toward' rounds toward zero
   - 'away' (default) rounds away from zero
   - 'up' rounds up
   - 'down' rounds down
   - 'even' rounds toward the even integer
   - 'odd' rounds toward the odd integer


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] polynomials over QQ

[D. S. McNeil]

 My understanding was that 'x' was the indeterminate
 of the ring of polynomials over QQ, i.e. the rationals.  So how come
 the polynomial
 has coefficients which are not rational?

Because the polynomial isn't living where you think it does anymore:

sage: R.x = QQ['x']
sage: R
Univariate Polynomial Ring in x over Rational Field
sage: parent(x)
Univariate Polynomial Ring in x over Rational Field
sage:
sage: f = (x-sqrt(2))*(x+sqrt(2))
sage: f
(x - sqrt(2))*(x + sqrt(2))
sage: parent(f)
Symbolic Ring
sage: R(f)
---
TypeError Traceback (most recent call last)
[...]
TypeError: unable to convert -sqrt(2) to a rational
sage: R(expand(f))
x^2 - 2

In this case, the introduction of the sqrt terms pushed the expression
out of R and into SR.  We can convert back, but only if the expression
is in a form that Sage can recognize as belonging to R.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: ValueError: setting an array element with a sequence. when plotting R data

[D. S. McNeil]

 Thus I wonder why this doesn't work

 def f(x): return r.dnorm(x,mean=100,sd=25)._sage_()
 plot(f(x)(x,90,110))

 while

 def f(x): return x^3
 plot(f(x),(x,0,2))

 works perfectly.

The problem is that the f(x) calls f with the symbolic argument x
_at the time you call the plot function_.

In the first case, r.dnorm doesn't know what to do with a symbolic x,
and barfs before it returns.  OTOH, in the second function, the f(x)
call simply returns x^3, so it's equivalent to typing plot(x^3, (x, 0,
2).

Both cases should work if you use

plot(f, (x,90,110))
plot(f, (x,0,2))


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] How can I tell if an algebraic number is rational?

[D. S. McNeil]

 PolynomialRing(ZZ, 'x')

This is only an aside, but I should probably warn that (unlike var,
say) this doesn't change x, so it might not do what you're thinking.

x is still an Expression, an element of the Symbolic Ring, and so f is
also an Expression.  You probably want to use something like

sage: R.x = PolynomialRing(ZZ)
sage: parent(x)
Univariate Polynomial Ring in x over Integer Ring
sage: f=x^3 + A*x +B
sage: parent(f)
Univariate Polynomial Ring in x over Integer Ring

or R.x = ZZ[] or something.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] weird inequality bug...

[D. S. McNeil]

The key here is understanding exactly what x[0] is: it's not a
rational.  If you run your code (after adding the line
set_random_seed(3) at the start to make sure we're working with the
same matrices), you see:

sage: minx, maxx, miny, maxy
(+Infinity, (3), +Infinity, (21/5))

and the odd parentheses should hint at the problem:

sage: type(maxx), type(maxy)
(type 'sage.modules.vector_rational_dense.Vector_rational_dense',
type 'sage.modules.vector_rational_dense.Vector_rational_dense')
sage: x[0]
(5/7)
sage: type(x[0])
type 'sage.modules.vector_rational_dense.Vector_rational_dense'

We're not comparing the rational 5/7 with infinity, we're comparing
the *vector* (5/7,) -- comma inserted for clarity -- with infinity.
Anyway, comparisons between things which shouldn't be compared tend to
give weird results if they work:

sage: x = 5/7
sage: x  infinity # good
True
sage: x  infinity # good
False
sage: x = vector([5/7])
sage: x  infinity
False
sage: x  infinity
True
sage: 5/7  vector([0.0])
True

If you replace x[0] by x[0][0], and x[1] by x[1][0], so that you're
comparing the entries and not the row vectors, it should do what you
expect.  (You could also use x = x.list() to coerce to a list and then
x[0] and x[1] will work, I guess, but that kind of hides what's going
on.)


Hope that helps,

Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Legend Options in 2D Plots? Title Options in 3D plots?

[D. S. McNeil]

 1 - On 2D plots, I have been unable to use the legend() attribute and
 the set_legend_options() attribute.

I don't think you have those methods.  I think you're using an old
version of sage (maybe = 4.5.2?) but reading documentation from a
later version, which occasionally causes trouble.  What does version()
report?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: How to evaluate an expression numerically

[D. S. McNeil]

 time[0].right().n()

 Not sure why time.rhs() doesn't work...

That I can explain.  time is a list (well, technically a Sequence) of
equations, and so you'd need to type time[0].rhs().  The individual
equations have right hand sides, but the list itself doesn't.

FWIW, I prefer using dicts rather than equations as the return type:

sage: var(A k c)
(A, k, c)
sage: phi(t) = A*exp(k*c*t)
sage: s = solve(phi==1,t,solution_dict=True)
sage: s
[{t: log(1/A)/(c*k)}]
sage: s[0][t]
log(1/A)/(c*k)

which handles the cases with more than one variable a little easier, I
think, but YMMV.


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Straightforward calculation going wrong?

[D. S. McNeil]

The range function is a Python one, and it returns Python ints.
Python ints have truncating division, so that 3/2 = 1, not 3/2. When
you type 3/2 at the Sage command, it's preparsed to be Sage Integers:

sage: 3/2
3/2
sage: preparse(3/2)
'Integer(3)/Integer(2)'
sage: int(3)/int(2)
1
sage: 3r/2r
1

(I don't know what the r is supposed to stand for but I always think
of it as raw, i.e. Python, not Sage.  I might even be right.)  So in
your code your C values are wrong in the second loop.   They work in
the first loop because B=-2 makes it a Sage Integer, and so the
results of divisions can become Rationals like you expect.


print 'Loop 2:'
for A in range(-3,-1):
   for B in range(-3,-1):
   Cint = -A*B/(A+B)
   CInt = (1)*-A*B/(A+B)
   print A, B,
   print 'Cint:', Cint, A*B+Cint*(A+B),
   print 'CInt:', CInt, A*B+CInt*(A+B)

gives:

-3 -3 Cint: 1 3 CInt: 3/2 0
-3 -2 Cint: 1 1 CInt: 6/5 0
-2 -3 Cint: 1 1 CInt: 6/5 0
-2 -2 Cint: 1 0 CInt: 1 0

The multiplication by the Sage Integer 1 makes the CInt expression a
Sage one, and so it works.  For these reasons, I tend to avoid using
range in Sage code entirely.  There are several
alternatives.  There's srange/xsrange=sxrange

sage: srange(-3, -1)
[-3, -2]
sage: sxrange(-3, -1)
generator object generic_xsrange at 0x10d57ae10
sage: list(sxrange(-3, -1))
[-3, -2]

(The s stands for Sage):

Or you can be explicit and use IntegerRange, which gives a more
informative string:

sage: IntegerRange(-3, -1)
{-3, -2}

I like the (a..b) and [a..b] syntaxes myself:

sage: for A in (-3..-2): print A, type(A)
:
-3 type 'sage.rings.integer.Integer'
-2 type 'sage.rings.integer.Integer'

but note that it includes the right hand limit.  I find this is more
useful in mathematical code than in pure programming, mostly because I
often find cases where I want (a..a) to include a, but your mileage
may vary.

Does that make sense?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support]

[D. S. McNeil]

 Emailing this instead of newpost because of the attachment.
 This worksheet works well ONCE when it starts unevaluated. Reediting the
 data followed by ActionEvaluate all induces TypeError: 'str' object
 is not callable. My way out is to Delete output save and quit and
 reopen. But I do not understand what happens.

Even without looking at the code, from the error message alone, I know
that there's a function name that you've rebound to a string, and when
you're running it the second time, it's got the string value instead
so the call fails.  So the only question is which function it is!

Looking at the code quickly reveals the problem: you've reassigned the
name n, but you also use it as an abbreviation for numerical_approx,
so that the second time it has the value Jeu 07h00, which is a
string, so n(mean((Sys[3*i],Sys[3*i+1],Sys[3*i+2])),digits=3) and
similar calls fail.

Workarounds (choose your favourite):

(1) add a reset() at the start of the worksheet, so everything returns
to normal, including n.  This doesn't solve the problem that after
you use n as a variable you can't use it for numerical_approx,
though.

(2) Use numerical_approx instead of n.

(3) Avoid n as a variable.

(4) Choose something else -- nn, or nap, or whatever -- and use it
instead (nn = numerical_approx, and then nn(3.2445,digits=2).

BTW, I see you're still writing

for i,j in enumerate(temps):
for k,l in enumerate(mesuresF):
for m,n in enumerate(heuresF):
if m==k:
if i==k:
print l[-1],  ,l[0],  ,l[1:4]

It's more Pythonic to write something like

for t,m,h in zip(temps, mesuresF, heuresF):
print m[-1],  ,m[0],  ,m[1:4]

After all, if there are a thousand data points in
temps/mesuresF/heuresF, then the above loop will take forever to
finish (1000^3=10^9)!

Flat is better than nested, we say, although in this particular loop
you only seem to be using one of the variables anyway, so for m in
mesuresF: looks like it should work.

For those of us who grew up with C or Fortran, it's hard to get out of
the habit of using loop indices everywhere, but it makes life much
better once you do. :^)

[For the record, I think encouraging people to use n as an alias for
numerical_approx was an unfortunate decision in the first place; it's
about a thousand times too common as a variable name to use as a
function name.]

Hope that helps!


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: memory error

[D. S. McNeil]

 TypeError: list indices must be integers, not
 FiniteField_ext_pariElement

 I see this problem on Mac OS X, so it's not specific to one type of system.  
 It may be due to a change in the way finite fields are handled as the size 
 grows, but I'm not familiar with that code.

 Can anyone hazard a guess?

I think there are two issues:

(1) _points_cache_sqrt, which is called if the order of the base ring
is = zech_log_bound=65536, requires the element to have an __index__
method.  This is what's giving the error message.  If I replace the
list with a dict and check for membership (or, I guess, the same
try/except KeyError approach as in _points_fast_sqrt) so __index__ is
unnecessary, everything seems to work.

The duplication in _points_fast_sqrt/_points_cache_sqrt suggests it
could do with a refactor.

.. but (2) Isn't this going to have an *extremely* large number of
points and take at least a month to compute?  IOW, even if we fix what
I think is the problem, I'm not sure if it's going to be very good at
doing what the OP wants it to.  I could easily be wrong, though, I'm
just naively extrapolating..


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] numpy divide by zero warnings are annoying

[D. S. McNeil]

 I'm not turning off warnings in numpy, though, since we use it under the hood 
 only
 here.

I'm confused.  I was going to recommend numpy.seterr(all='ignore')
before I read this, maybe wrapping plot to restore the original state
after the call..  but now I'm not sure what kind of solution you want.
 :^)

In this case, the problem is being caused by the default parameters in
plot_slope_field.  If you override headlength=0 with some small
number, there's no problem.  IOW,

P=plot_slope_field(g,(x,3,4),(y,-1,1),headlength=1e-8)

works for me.  FYI, it's the following few lines in Quiver._h_arrows at fault:

minsh = self.minshaft * self.headlength
[]
shrink = length/minsh
X0 = shrink * X0[np.newaxis,:]
Y0 = shrink * Y0[np.newaxis,:]

Probably we should change the defaults and/or (if it's not done
already) ask our matplotlib friends to special-case 0 for no
arrowheads.

Does that help?


Doug

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Using scipy special functions

[D. S. McNeil]

 I'm completely unable to get the scipy special functions module to
 work. In addition, it seems to cause chaos on my system once imported


 sage: import scipy
 sage: from scipy.special import *
 sage: scipy.special.lpn(1,1)

I'd avoid your second line in Sage, which pulls everything in
scipy.special into the local namespace.  The *-import means that a lot
of Sage functions are replaced by (num/sci)py identically-named
functions, which will cause no end of trouble.  (E.g. the standard
Sage symbolic cos is replaced by the numpy cos.)  I'm guessing that's
the chaos you're referring to.

I'd recommend using

  import scipy.special

especially if you're going to be writing scipy.special.lpn(1,1) anyway.

As for the problem itself, this is a known issue, due to an
unfortunate design decision in the numpy/scipy isscalar function.  It
which doesn't duck type, and so objects to types it doesn't know such
as those used by Sage (or even Python-native types not in its list).

http://trac.sagemath.org/sage_trac/ticket/10928

Lots of the scipy.special functions _do_ work, but any that use
isscalar will probably fail.  A workaround is to make sure that what
lpn is passed is a Python object:

sage: scipy.special.lpn(int(1), float(2))
(array([ 1.,  2.]), array([ 0.,  1.]))
sage: scipy.special.lpn(1r, 2r)
(array([ 1.,  2.]), array([ 0.,  1.]))

    597     if not (isscalar(n) and isscalar(z)):
 -- 598         raise ValueError, arguments must be scalars.
    599     if (n!= floor(n)) or (n0):
    600         raise ValueError, n must be a non-negative integer.

 ValueError: arguments must be scalars.
 sage: n
 function numerical_approx at 0x2dc1668

Just for the record, this n isn't the n in the code fragment, which
is the first argument to lpn.


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Qhull package installation fail

[D. S. McNeil]

 I tried installing the Qhull package which is an optional package on
 the list at http://www.sagemath.org/packages/optional/, but i got the
 error below.  It appears that /usr/include/float.h can't be found, but
 i *do* have that file.

I'm not sure whether this is a Sage-related issue.  Could be revealing
a preexisting config problem on your system in your gcc install.  The
relevant line of float.h appears to be

#include_next float.h

which tries to find the *next* float.h in the include path after the
current one, and it's that second float.h that it's failing to find.
But I'm not sure why gcc is looking in /usr/include/float.h.  I think
it should pick up your equivalent of
/usr/lib/gcc/i686-apple-darwin10/4.2.1/include/float.h first instead,
and never make it to /usr/include/float.h in the first place.  It
feels kind of like gcc doesn't know where it lives.

Does the code

#include float.h
int main(){return 0;}

compile for you?  [gcc -E and gcc -v output if possible.]


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] edge labels in trees

[D. S. McNeil]

 I want to show a tree  with edge labels. I tried out

 G1 = Graph({1:{5:0},2:{5:1},3:{6:1},4:{7:0},5:{6:0},6:{7:1}})
 show(G1,edge_labels=true)
 show(G1,layout=tree,edge_labels=true)

 In the first graphic the labels are positioned on the edges but in the
 second one they are not. Why is this? How can I get it right?

I think this is related to
http://trac.sagemath.org/sage_trac/ticket/10124, or at least my
patched version seems to work.  The basic problem is that
graph_plot.py has problems with integer positions because integer
division is truncating in Python 2 and it's easy to forget to handle
the integer case when your inputs are typically floats.  (This problem
goes away in Python 3.)

You could either apply the posted patch, or -- as a hack which doesn't
require any rebuilding -- you could type the following before you plot
G1:


# trac #10124 partial monkeypatch
def set_pos(self):
self._pos = self._graph.layout(**self._options)
self._pos = dict((k,(float(v[0]), float(v[1]))) for k,v in
self._pos.iteritems())

import sage.graphs.graph_plot
sage.graphs.graph_plot.GraphPlot.set_pos = set_pos


which should work around the problem in the meantime, unless this
coffee hasn't kicked in yet..


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] timing out a function or command

[D. S. McNeil]

On Tue, Apr 5, 2011 at 1:47 PM, tvn nguyenthanh...@gmail.com wrote:
 Is there a timeout mechanism in Sage or Python that allows me to kill or
 raise an exception on a command that exceeds some time threshold ?

You probably want alarm, which raises a KeyboardInterrupt.


Doug

-- 
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] using random seed for sample() function

[D. S. McNeil]

On Tue, Mar 29, 2011 at 4:00 AM, tvn nguyenthanh...@gmail.com wrote:
 I'd like to be able to regenerate  samples by feeding a seed value to
 random.seed()  ,  but it seems sample() doesn't use this random seed.  Is
 there a way to do what I want ?

help(sage.misc.randstate) explains a lot of the gory details.

set_random_seed() controls the _Sage_ random seed, but I think if
you're replacing the random() function with the Python random module
(to use random.random() instead of just random()) you're going to have
to set both:

sage: reset()
sage: import random
sage:
sage: # try just set_random_seed..
sage: set_random_seed(1)
sage: random.random()
0.84743373693723267
sage: sample(range(5),2)
[4, 2]
sage:
sage: set_random_seed(1)
sage: random.random()
0.76377461897661403
sage: # that didn't work..
sage: sample(range(5),2)
[4, 2]
sage: # but this did!
sage:
sage: # so we have to use both..
sage: set_random_seed(1) # for Sage
sage: random.seed(1) # for Python random module
sage: random.random()
0.13436424411240122
sage: sample(range(5),2)
[4, 2]
sage:
sage: set_random_seed(1) # for Sage
sage: random.seed(1) # for Python random module
sage: random.random()
0.13436424411240122
sage: sample(range(5),2)
[4, 2]


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Units Package Simplification in Exponentiation/Powers

[D. S. McNeil]

 sage: (a^3)^(1/3)
 5*(meter^3)^(1/3)

 does not produce the expected units of meter.  Is there a means to
 force further simplification?

You could try

sage: m = units.length.meter
sage: assume(m  0)
sage: (5*m)^3
125*meter^3
sage: ((5*m)^3)^(1/3)
5*(meter^3)^(1/3)
sage: simplify(((5*m)^3)^(1/3))
5*meter



Doug, who likes his Yankee friends but insists that 'meter' is a
device which measures ;-)

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Additionnal constraints are never enough

[D. S. McNeil]

 integrate(integrate(1/16,y,-2,(z+2*x+4)/(x+2)),x,(-z-8)/4,2)

Seems to work for me:

--
| Sage Version 4.6.1, Release Date: 2011-01-11   |
| Type notebook() for the GUI, and license() for information.|
--
sage: var(x y z)
(x, y, z)
sage: integrate(integrate(1/16,y,-2,(z+2*x+4)/(x+2)),x,(-z-8)/4,2)
[...]
TypeError: Computation failed since Maxima requested additional
constraints (try the command 'assume(z+160)' before integral or limit
evaluation, for example):
Is  z+16  positive, negative, or zero?
sage: assume(z+160)
sage: integrate(integrate(1/16,y,-2,(z+2*x+4)/(x+2)),x,(-z-8)/4,2)
[...]
TypeError: Computation failed since Maxima requested additional
constraints (try the command 'assume(z0)' before integral or limit
evaluation, for example):
Is  z  positive, negative, or zero?
sage: assume(z0)
sage: integrate(integrate(1/16,y,-2,(z+2*x+4)/(x+2)),x,(-z-8)/4,2)
-1/16*z*log(-1/4*z) + 1/16*z*log(4) + 1/16*z + 1

What version of Sage are you running?  (Note that 4.6.2 has just been released.)


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Additionnal constraints are never enough

[D. S. McNeil]

 But, after a reset() command, it does not work anymore.

Confirmed.  Well, that ain't proper.

Could you try reset(); forget(); instead of reset()?  Right now
assumptions survive a reset in a broken state (try assumptions()
before and after the reset and you'll see them still listed even if
the variable isn't defined any longer) and I wonder if that's causing
problems.  I have a positively-reviewed patch --
http://trac.sagemath.org/sage_trac/ticket/10855 -- pending for 4.7
which fixes this.  I think there are other ways to get the maxima
interface into an unhappy state, but this works too.

--
| Sage Version 4.6.1, Release Date: 2011-01-11   |
| Type notebook() for the GUI, and license() for information.|
--

sage: var(x y z)
(x, y, z)
sage: assume(z+160, z0)
sage: integrate(integrate(1/16,y,-2,(z+2*x+4)/(x+2)),x,(-z-8)/4,2)
-1/16*z*log(-1/4*z) + 1/16*z*log(4) + 1/16*z + 1

sage: reset()
sage: z
[...]
NameError: name 'z' is not defined
sage: assumptions()
[z + 16  0, z  0]   # - bug!, should be []

sage: var(x y z)
(x, y, z)
sage: assume(z+160, z0)
sage: integrate(integrate(1/16,y,-2,(z+2*x+4)/(x+2)),x,(-z-8)/4,2)
ERROR: An unexpected error occurred while tokenizing input
[...]
TypeError: Computation failed since Maxima requested additional
constraints (try the command 'assume(z+160)' before integral or limit
evaluation, for example):
Is  z+16  positive, negative, or zero?

sage: forget()
sage: assume(z+160, z0)
sage: integrate(integrate(1/16,y,-2,(z+2*x+4)/(x+2)),x,(-z-8)/4,2)
-1/16*z*log(-1/4*z) + 1/16*z*log(4) + 1/16*z + 1


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: modifying a plot after creation

[D. S. McNeil]

 Note that you need to use randomize=False in plot if you're doing
 this, because otherwise the plot positions are different each time
 (!), so caching the function is useless.  I don't understand the
 motivation for having that be the default behaviour.

 Better than having a set algorithm that will perform poorly in edge
 cases (e.g. periodic functions).  The adaptive algorithm could
 certainly be improved - it nearly always produces more points than you
 ever need, and produces monstrosities like plot(tan(x),-10,10) - but
 overall this has served us well.

Adaptive != non-reproducible, though.

I can understand why having an adaptive approach which uses random
numbers is more robust than dx=(xmax-xmin)/100 etc.  What I don't
understand is why the default behaviour isn't to use a set seed so
that the result of plotting calls would be reproducible.

If there's a real difference between using one seed and the next, in
that sometimes it produces a good plot and sometimes it doesn't, then
aren't you going to need to control the seed anyway to get the good
plot back?  And if there isn't a real difference, then why not simply
pick one and get all the benefits of the adaptive approach while
preserving reproducibility?

(Having re-read the docs it looks like I should have added
adaptive_recursion=0, not just randomize=False.)

To avoid possibly causing problems for tough functions, then, I guess
setting the random seed is better than using the randomize flag:

# or a python function:
def f(x):
print 'called with',x
return x**2

g = CachedFunction(f)
set_random_seed(3)
plot(g, 0, 10) # this is very noisy because of the print
set_random_seed(3)
plot(g, 0, 10) # this is quiet


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] sage upgrade on MacOs, Japanese envirinment.

[D. S. McNeil]

On Wed, Mar 2, 2011 at 4:55 PM, Thierry Dumont
tdum...@math.univ-lyon1.fr wrote:
 One of my colleagues uses a Macintosh with a Japanese environment.
 He cannot upgrade sage (4.6.1).Here is a transcription of what happens.
 Any idea?

I used to hit this problem myself (not with Sage, but with other
software) because of http://bugs.python.org/issue1276.

In this case I think entering

export HGENCODING=UTF-8

from the command line might fix the problem (at least it did for me,
if I successfully mocked up the circumstances).


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] sage upgrade on MacOs, Japanese envirinment.

[D. S. McNeil]

 export HGENCODING=UTF-8

(Of course not everyone uses bash, so modify as appropriate.)


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Counting lattice points inside convex polyhedra

[D. S. McNeil]

 I was wondering if sage implements any algorithm for counting the number of
 points with integer coordinates inside polyhedra with rational coordinates.
  even such an algorithm for polygons would be useful for me.

Have a look at the integral_points method of Polyhedron objects, which
might do what you're looking for:

sage: P = Polyhedron([[1/2, 1/2], [1/2, 7/2], [7/2, 1/2], [7/2, 7/2]])
sage: P
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices.

sage: P.integral_points()
[(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3)]
sage: len(P.integral_points())
9

Note that this constructs the integral points rather than counts them,
so if there are lots of points in the enveloping lattice polytope it
won't be that fast, but maybe it'll work for your purposes.


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Problem integrating piecewise linear function (breaks with constant piece)

[D. S. McNeil]

On Thu, Feb 24, 2011 at 1:52 PM, Ronald L. Rivest rivest@gmail.com wrote:
 Is posting this bug here on this list sufficient to report it so that it
 will (eventually) get fixed, or is there some other process for doing so
 that needs to be done next?

I think in general you open a ticket on trac.
http://trac.sagemath.org/sage_trac/
I find http://trac.sagemath.org/sage_trac/report/33 [recent activity]
the most useful
front page.

In this case I think I see what the problem is, so I've opened a
ticket and put up a patch,
#10841.


Doug

-- 
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Numeric multiple integral of expression not integrable symbolically ?

[D. S. McNeil]

On Thu, Feb 24, 2011 at 1:14 PM, Ronald L. Rivest rivest@gmail.com wrote:

 I want a numeric integration, but
numerical_integral(numerical_integral(f,0,x),0,1)
 doesn't work, since f takes two parameters, not one.

IIUY, you could nest it, and do

sage: numerical_integral(lambda x: numerical_integral(lambda y:
(y+1)^(y+1)^x,(0,x))[0],(0,1))
(0.76066264344246259, 8.4450518072205787e-15)

or

sage: import mpmath
sage: f = mpmath.quad(lambda x: mpmath.quad(lambda y: (y+1)^ (y+1)^
x,(0,x)),(0,1))
sage: f
mpf('0.76066264344246248')

I like mpmath because I more often want the increased accuracy than I
care about the decrease in speed.

sage: mpmath.mp.dps = 50
sage: f = mpmath.quad(lambda x: mpmath.quad(lambda y: (y+1)^ (y+1)^
x,(0,x)),(0,1))
sage: f
mpf('0.76066264344246243183604278581374359729636608721856373')

But since flat is better than nested, I think something like

numerical_integral(f, (x,0,3), (y,0,x^2))

should Just Work(tm).


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Re: Isolating real roots of exact univariate polynomial

[D. S. McNeil]

A somewhat simpler test case, which I think preserves the qualitative issue:

sage: from sage.rings.polynomial.real_roots import real_roots
sage:
sage: x = polygen(QQ)
sage: f = 2503841067*x^13 - 15465014877*x^12 + 37514382885*x^11 -
44333754994*x^10 + 24138665092*x^9 - 2059014842*x^8 - 3197810701*x^7 +
803983752*x^6 + 123767204*x^5 - 26596986*x^4 - 2327140*x^3 + 75923*x^2
+ 7174*x + 102
sage: len(real_roots(f)), len(real_roots(f,strategy='warp'))
(11, 13)
sage:
sage: [map(float,r[0]) for r in real_roots(f)]
[[-0.28173828138369572, -0.23901367199914603], [-0.19628906261459633,
-0.15356445323004664], [-0.078796386807084673, -0.073455810634015961],
[-0.036071777422534979, -0.030731201249466267],
[-0.025390625076397555, 0.017333984308152139], [0.060058593692701834,
0.10278320307725153], [0.14550781246180122, 0.23095703123090061],
[0.57275390625, 0.615478515625], [0.658203125, 0.76618289947509766],
[1.375, 1.4375], [1.5, 2.0]]
sage: [map(float,r[0]) for r in real_roots(f,strategy='warp')]
[[-0.1, -0.23076923076923075], [-0.18518518518518517,
-0.14285714285714285], [-0.10344827586206896, -0.06],
[-0.049180327868852458, -0.032258064516129031],
[-0.023997, -0.015873015873015872], [0.032258064516129031,
0.06], [0.14285714285714285, 0.23076923076923075],
[0.59998, 0.64102564102564097], [0.72972972972972971,
0.75342465753424648], [0.95373241116145957, 0.95396541443053062],
[0.95419847328244267, 0.96923076923076923], [1.2857142857142856,
1.4615384615384615], [1.6665, 3.0]]


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Numpy Troubles

[D. S. McNeil]

I think it's neither Sage nor numpy that's at fault, it's a weird
interaction (Sage 4.6.1):

sage: import numpy
sage: numpy.binary_repr(17)
''
sage: numpy.binary_repr(int(17))
'10001'

and I think it's rooted in this fact:

sage: hex(17)
'11'
sage: hex(int(17))
'0x11'

That is, Sage capital-I Integers deliberately don't have the '0x'
prepended, and the numpy.binary_repr routine has lines

ostr = hex(num)
[..]
bin = ''.join([_lkup[ch] for ch in ostr[2:]])

which fail if the prefix isn't there.


Doug

-- 
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] series Bug

[D. S. McNeil]

 Sorry for the large example, but smaller examples seem to work.

After some reductions, I think we can find a slightly simpler failure,
of a somewhat suspicious size (4.6.1, OS X 10.6):

q = var('q')
for i in [1..33]:
ix = 2**i
f=(q+1)/(q^(ix) + 1)
res = (f.series(q==0,2), f.subs(q=0), f.taylor(q, 0, 2))
print i, ix, res

1 2 (1 + 1*q + Order(q^2), 1, -q^2 + q + 1)
2 4 (1 + 1*q + Order(q^2), 1, q + 1)
3 8 (1 + 1*q + Order(q^2), 1, q + 1)
4 16 (1 + 1*q + Order(q^2), 1, q + 1)
[...]
28 268435456 (1 + 1*q + Order(q^2), 1, q + 1)
29 536870912 (1 + 1*q + Order(q^2), 1, q + 1)
30 1073741824 (1 + 1*q + Order(q^2), 1, q + 1)
31 2147483648 (Order(q^(-2147483646)), 1, q + 1)
32 4294967296 (1 + 1*q + Order(q^2), 1, q + 1)
33 8589934592 (1 + 1*q + Order(q^2), 1, q + 1)

The case with q^(3^i) is even funnier:

18 387420489 (1 + 1*q + Order(q^2), 1, q + 1)
19 1162261467 (1 + 1*q + Order(q^2), 1, q + 1)
20 3486784401 (Order(q^2), 1, q + 1)
21 10460353203 (1 + 1*q + Order(q^2), 1, q + 1)
22 31381059609 (1 + 1*q + Order(q^2), 1, q + 1)
23 94143178827 (Order(q^2), 1, q + 1)
24 282429536481 (Order(q^2), 1, q + 1)
25 847288609443 (1 + 1*q + Order(q^2), 1, q + 1)
26 2541865828329 (Order(q^2), 1, q + 1)
27 7625597484987 (1 + 1*q + Order(q^2), 1, q + 1)
28 22876792454961 (1 + 1*q + Order(q^2), 1, q + 1)
29 68630377364883 (1 + 1*q + Order(q^2), 1, q + 1)
30 205891132094649 (Order(q^2), 1, q + 1)
31 617673396283947 (1 + 1*q + Order(q^2), 1, q + 1)
32 1853020188851841 (Order(q^2), 1, q + 1)
33 555906056623 (Order(q^2), 1, q + 1)

Obvious suspect is overflow due to use of a C int somewhere along the
path (suddenly negative around 2^31; works for small numbers; is
better-behaved for powers of 2; can generate an actual OverflowError
at large enough exponents.)


Doug

-- 
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] plotting box function

[D. S. McNeil]

 hello,
 why is the below code plotting a flat function rather than a box one?

There are two things going on.  First, in the line

plot(box(x,1),(x,-3,3))

box(x,1) is actually being evaluated when the line is executed, and
not thereafter.  IOW you're computing box(x, 1), which is 0, so the
above is equivalent to

plot(0, (x, -3, 3))

You might want to stick a print statement inside the box function
(print called!) to convince yourself this is true.

Second, box(x,1) = 0 because the condition abs(x)  1 is False for a
variable x, and so the else is executed.  Note that False here
translates as I can't prove that it's True: if instead of the else
you'd written abs(x) = 1 ,that'd be False too, neither path would
get executed, and so the result would be None (what Python returns
when there's no explicit return statement.)  OTOH, if you call box(x,
infinity), you get 1.

There are a few ways around this.  Probably the most general-purpose
solution (which works even when some Sage-specific tricks don't) is to
delay the execution of the box function by writing a lambda-function
wrapper:

plot(lambda x: box(x,1), (x, -3, 3))

which is a short way to avoid having to write a new function def
box1(x): return box(x, 1) and then calling plot(box1, (x, -3, 3)).


Doug

-- 
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Re: [sage-support] Plotting: axes with ticks but with no labels

[D. S. McNeil]

On Fri, Jan 28, 2011 at 9:08 PM, Jeff wrote:
 I would like to be able to plot a function, e.g. plot(sin), that has
 axes and ticks on the axes but that does not have labels for the
 ticks. I understand that I might be able to do this using a ticker
 formatter, perhaps also, by directly using matplotlib, but I do not
 know exactly how to go about doing this.

There may be a simpler way, but:

import matplotlib

p = plot(sin)
p.show(tick_formatter=(matplotlib.ticker.NullFormatter(),
   matplotlib.ticker.NullFormatter()))

worked for me.  The repetition is to make sure that both x and y tick
labels are turned off.


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

[sage-support] Re: Precision in (Real) Vector Space / Working with reflection groups

[D. S. McNeil]

 (-6.66133814775094e-16, -1.00)
 ---???

 Why isn't the first coordinate 0??? In theory it should be...

Such things are basically inevitable with floating point whenever any
number involved can't be perfectly represented-- and in this case
sqrt(2) is irrational to boot.  Sometimes you luck out and sometimes
you don't.  help(RR) gives a brief warning.  Here I think you can
avoid the problem by working not in RR (2^(-53) ~ 1e-16) but the
symbolic ring SR.

sage: V=VectorSpace(RR,2,inner_product_matrix=[[1, -sqrt(2)/2], [-
sqrt(2)/2, 1]])
sage: V.basis()
[
(1.00, 0.000),
(0.000, 1.00)
]
sage: V=VectorSpace(SR,2,inner_product_matrix=[[1, -sqrt(2)/2], [-
sqrt(2)/2, 1]])
sage: a,b = V.basis()
sage: a,b
((1, 0), (0, 1))
sage: r1=Sigma(a,b);r1
(sqrt(2), 1)
sage: r2=Sigma(r1,a);r2
(-1, -sqrt(2))
sage: r3=Sigma(r2,r1);r3
(0, -1)

If you really want to work in floating point you're going to have to
loosen your definition of equality -- see, e.g.,
http://floating-point-gui.de/errors/comparison/ -- possibly doing so
rigorously with interval arithmetic.


Doug

--
Department of Earth Sciences
University of Hong Kong

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org