memory leak?

[opti]

Hello,

I am using sage 3.4
I have written a simple sage script to evalute the cardinal of Brent-
Suyama elliptic curves s=11
for primes < 2^27.

Memory usage increases till program (or machine) crashes.

The script is the following:

p=30

def FindGroupOrder(p,s):
K = GF(p)
v = K(4*s)
u = K(s**2-5)
x = u**3
b = 4*x*v
a = (v-u)**3*(3*u+v)
A = a/b-2
x = x/v**3
b = x**3 + A*x**2 + x
E = EllipticCurve(K,[0,b*A,0,b**2,0])
return factor(E.cardinality())

while p<134217689:
p=next_prime(p)
g=FindGroupOrder(p,11)
print g


Any idea why this happens?

Regards,

[mabshoff]

On Apr 29, 11:11 am, opti <iram.che...@gmail.com> wrote:
> Hello,

Hi,

I would guess it is a leak :)

Running it up to p=966233 consumed about 2 GB - I will have a closer
look. IIR there are some open memory leaks in that area of the code,
but I should know more in a couple hours.

> Regards,

Cheers,

Michael

[opti]

Michael,
Thanks for your reply.
I did the exact same thing under magma, and it consumes 14 meg.
Have you found anything new?
Regards,

Iram

On Apr 29, 9:07 pm, mabshoff <Michael.Absh...@mathematik.uni-

[mabshoff]

On Apr 30, 9:15 am, opti <iram.che...@gmail.com> wrote:
> Michael,

Hi Iram,

> Thanks for your reply.
> I did the exact same thing under magma, and it consumes 14 meg.

Hehe, you know how to get us going :)

> Have you found anything new?

I did find the cause, but I have not found a fix yet. What seems to
happens is that for each p a multivariate polynomial ring in
libSingular is leaked - all the gory details are at

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

> Regards,
>
> Iram

Cheers,

Michael

[John Cremona]

The actual cardinality is computed using the pari library (since the
field is a prime field and p<10**18). I the leak is in libsingular it
must come from constructing the curve, not from computing the
cardinality. So it would be worth testing the loop with a new
function which constructs the curve and returns nothing.

John

On Apr 30, 10:47 pm, mabshoff <Michael.Absh...@mathematik.uni-

[mabshoff]

On May 1, 2:02 am, John Cremona <John.Crem...@gmail.com> wrote:

Hi John,

> The actual cardinality is computed using the pari library (since the
> field is a prime field and p<10**18).  I the leak is in libsingular it
> must come from constructing the curve, not from computing the
> cardinality.  So it would be worth testing the loop with a new
> function which constructs the curve and returns nothing.

Ok, this seems to indicate for me that maybe Coercion could be
involved here in some shape or form since it seems that for every
prime we iterate over we create/leak a new mv ring in libSingular. And
that cannot be a good thing :)

No that I have a clue how to debug this, so I will poke RobertWB - I
am sure he has nothing better to do that to dive into this problem and
will be happy to debug this ;)

> John

Cheers,

Michael

[Robert Bradshaw]

Yep, I bet coercion is leaking rings. I'll look at this during Sage
days, if not sooner.

- Robert

[mabshoff]

On May 1, 12:12 pm, Robert Bradshaw <rober...@math.washington.edu>
wrote:

<SNIP>

> > No that I have a clue how to debug this, so I will poke RobertWB - I
> > am sure he has nothing better to do that to dive into this problem and
> > will be happy to debug this ;)
>
> Yep, I bet coercion is leaking rings. I'll look at this during Sage  
> days, if not sooner.

Thanks, I ought to look over your shoulder to get an idea on how to
attack those leaks in the future :)

> - Robert

Cheers,

Michael

[John Cremona]

When an elliptic curve is created the code in the __init__ function in
ell_generic.py (lines 164-5) do cause a multivariate polynomial ring
to be created. In this case it's a new ring each time as the base
field is always a new field.

The reason this is done is that the elliptic curve class derives from
the plane curve class, which is created from a homogeneous polynomial
in three variables. However at least 95% (at a guess) of elliptic
curve functionality does not depend at all on the base class (one
exception is checking that points lie on the curve, and even that was
not the case until the recent patch by Alex Ghitza).

So, when you create a large number of curves over different fields you
do also create all these rings. However, in the example in the
original post, these elliptic curves do not remain in scope, so I
don't know why the rings are not garbage-collected away.

John

On May 1, 10:00 pm, mabshoff <Michael.Absh...@mathematik.uni-

[simon...@uni-jena.de]

Dear John,

On 2 Mai, 18:39, John Cremona <John.Crem...@gmail.com> wrote:
> When an elliptic curve is created the code in the __init__ function in
> ell_generic.py (lines 164-5) do cause a multivariate polynomial ring
> to be created.  In this case it's a new ring each time as the base
> field is always a new field.

How is the multivariate polynomial ring created?

I remember that recently Martin Albrecht was fixing some bug: Here, a
multivariate polynomial ring was created without using the
PolynomialRing constructor, violating uniqueness of parents.

But if, in your case, you create a polynomial ring, work with it, and
then want it to be garbage collected, then perhaps it would be a good
idea to *not* use the PolynomialRing constructor: If the resulting
ring is cached (-> uniqueness of Parent Structures) then it would not
be garbage collected, or would it?

Best regards,
Simon

[John Cremona]

On May 2, 5:55 pm, simon.k...@uni-jena.de wrote:
> Dear John,
>
> On 2 Mai, 18:39, John Cremona <John.Crem...@gmail.com> wrote:
>
> > When an elliptic curve is created the code in the __init__ function in
> > ell_generic.py (lines 164-5) do cause a multivariate polynomial ring
> > to be created.  In this case it's a new ring each time as the base
> > field is always a new field.
>
> How is the multivariate polynomial ring created?

First the projective plane is constructed, effectively by
PP=ProjectiveSpace(2,K), and then the ring by PP.coordinate_ring(),
which calls
PolynomialRing(self.base_ring(),
self.variable_names(),
self.dimension_relative()+1)
so in effect we call PolynomialRing(K).

John

[mabshoff]

On May 2, 10:40 am, John Cremona <John.Crem...@gmail.com> wrote:
> On May 2, 5:55 pm, simon.k...@uni-jena.de wrote:
>
> > Dear John,
>
> > On 2 Mai, 18:39, John Cremona <John.Crem...@gmail.com> wrote:
>
> > > When an elliptic curve is created the code in the __init__ function in
> > > ell_generic.py (lines 164-5) do cause a multivariate polynomial ring
> > > to be created.  In this case it's a new ring each time as the base
> > > field is always a new field.
>
> > How is the multivariate polynomial ring created?
>
> First the projective plane is constructed, effectively by
> PP=ProjectiveSpace(2,K), and then the ring by PP.coordinate_ring(),
> which calls
>  PolynomialRing(self.base_ring(),
>                                self.variable_names(),
> self.dimension_relative()+1)
> so in effect we call PolynomialRing(K).
>
> John

Well, is there any way to avoid this? It seems like a bad thing to
waste a couple GB on rings that in the end aren't used any more once
the curve E is delete.

Cheers,

Michael

[Alex Ghitza]

>> > > When an elliptic curve is created the code in the __init__ function in
>> > > ell_generic.py (lines 164-5) do cause a multivariate polynomial ring
>> > > to be created.  In this case it's a new ring each time as the base
>> > > field is always a new field.
>>
>> > How is the multivariate polynomial ring created?
>>
>> First the projective plane is constructed, effectively by
>> PP=ProjectiveSpace(2,K), and then the ring by PP.coordinate_ring(),
>> which calls
>>  PolynomialRing(self.base_ring(),
>>                                self.variable_names(),
>> self.dimension_relative()+1)
>> so in effect we call PolynomialRing(K).
>>
>> John
>
> Well, is there any way to avoid this? It seems like a bad thing to
> waste a couple GB on rings that in the end aren't used any more once
> the curve E is delete.
>

I can confirm (if there was any doubt) that one can try

{{{
for p in prime_range(upper):
R.<x, y, z> = PolynomialRing(GF(p), 3)
}}}

and watch it eat 800MB of RAM.

I can see good reasons why we want the current behaviour of
"uniqueness of parents", and what we just saw here is a good reason
why we might *not* want it. How can we reconcile this? Can we make
it easy to turn off "uniqueness of parents" in some situations?

Best,
Alex

--
Alex Ghitza -- Lecturer in Mathematics -- The University of Melbourne
-- Australia -- http://www.ms.unimelb.edu.au/~aghitza/

[Robert Bradshaw]

We can have uniqueness and collection--it's called weakrefs. I
suspect part of the problem is that the coercion model doesn't use
them enough (there's places I meant to go back to and add them in).

- Robert

[simon...@uni-jena.de]

Hi!

On 3 Mai, 05:42, Alex Ghitza <aghi...@gmail.com> wrote:
...

> I can see good reasons why we want the current behaviour of
> "uniqueness of parents", and what we just saw here is a good reason
> why we might *not* want it.  How can we reconcile this?  Can we make
> it easy to turn off "uniqueness of parents" in some situations?

In a dirty way, certainly: The polynomial rings constructed with
PolynomialRing are cached at
sage.rings.polynomial.polynomial_ring_constructor._cache

So, just delete appropriate items in that cache.

Robert mentioned weak references. Does it mean the following?
- When creating a reference to a polynomial ring in that cache, then
tell the garbage collector that this reference should perhaps delay
the deletion of the polynomial ring, but it should not prevent
deletion.
- When there is any other reference to the polynomial ring, then this
counts as a proper reference, preventing the garbage collection.

Hence, when you do
sage: R=QQ['x']
then there is a weak reference in cache and a proper reference ``R``
to QQ['x']. So, the ring won't be garbage collected. When you continue
sage: R=GF(3)['x']
then only the weak reference points to QQ['x']. So, the ring may be
garbage collected.

I think in this way one would have unique parents *and* would fix the
memory leak, but I am afraid I don't know enough about weak references
to implement it.

Best regards,
Simon

[Robert Bradshaw]

Yes, that is exactly how it works. We should be using a
weakref.WeakValueDictionary (see http://docs.python.org/library/
weakref.html).

Perhaps changing line 39 of
sage.rings.polynomial.polynomial_ring_constructor from

_cache = {}

to

_cache = weakref.WeakValueDictionary

would be enough. For some reason, weakrefs used to be used, but
that's been disabled.

http://hg.sagemath.org/sage-main/diff/b4f33dc4908b/sage/rings/
polynomial/polynomial_ring_constructor.py

I think the underlying issue has been fixed. Note, however, that
there may also be caching elsewhere (e.g. I suspect a spot or two in
the coercion model).

- Robert

[simon...@uni-jena.de]

Hi!

On 3 Mai, 08:48, simon.k...@uni-jena.de wrote:
,,,

> I think in this way one would have unique parents *and* would fix the
> memory leak, but I am afraid I don't know enough about weak references
> to implement it.

Thanks to google, I found how it works:
Replace the line
_cache = {}
in sage.rings.polynomial.polynomial_ring_constructor by
from weakref import WeakValueDictionary
_cache = WeakValueDictionary({})


One can then do
sage: for p in primes(2,1000000):
....: R.<x,y,z> = GF(p)[]
sage: get_memory_usage()
778.35546875

Hmm. Are these the 800 MB of RAM that Alex was talking about?

Anyway. There is a patch awaiting review at # 5970


Cheers,
Simon

[simon...@uni-jena.de]

PS:

> One can then do
>   sage: for p in primes(2,1000000):
>   ....:     R.<x,y,z> = GF(p)[]
>   sage: get_memory_usage()
>   778.35546875

And:
sage: len(sage.rings.polynomial.polynomial_ring_constructor._cache)
34

Hence, indeed, some rings were garbage collected.

Cheers,
Simon

[mabshoff]

On May 3, 12:28 am, simon.k...@uni-jena.de wrote:

Hi,

Is there any reason you opened a new ticket and did not use #5949 as
mentioned above?

I am not sure this fixes the problem mentioned above, but I am still
testing. It might reduce the memory used already, but I won't know for
a while.

> Cheers,
>       Simon

As is this patch breaks badly:

sage -t -long devel/sage/sage/rings/number_field/
number_field.py # Segfault
sage -t -long devel/sage/sage/rings/tests.py # Segfault
sage -t -long devel/sage/sage/rings/number_field/
number_field_rel.py # Segfault
sage -t -long devel/sage/sage/rings/number_field/
number_field_element.pyx # Segfault
sage -t -long devel/sage/sage/rings/residue_field.pyx #
Segfault
sage -t -long devel/sage/sage/rings/number_field/
number_field_ideal_rel.py # Segfault
sage -t -long devel/sage/sage/rings/number_field/morphism.py #
Segfault
sage -t -long devel/sage/sage/rings/polynomial/
polynomial_singular_interface.py # Segfault
sage -t -long devel/sage/sage/rings/number_field/unit_group.py
# Segfault
sage -t -long devel/sage/sage/rings/number_field/
small_primes_of_degree_one.py # Segfault
sage -t -long devel/sage/doc/en/bordeaux_2008/nf_orders.rst #
Segfault
sage -t -long devel/sage/sage/rings/number_field/maps.py #
Segfault
sage -t -long devel/sage/sage/schemes/generic/affine_space.py
# Segfault

[Robert Bradshaw]

On May 3, 2009, at 1:04 AM, mabshoff wrote:

> On May 3, 12:28 am, simon.k...@uni-jena.de wrote:
>
> Hi,
>
> Is there any reason you opened a new ticket and did not use #5949 as
> mentioned above?
>
>> PS:
>>
>>> One can then do
>>> sage: for p in primes(2,1000000):
>>> ....: R.<x,y,z> = GF(p)[]
>>> sage: get_memory_usage()
>>> 778.35546875
>>
>> And:
>> sage: len(sage.rings.polynomial.polynomial_ring_constructor._cache)
>> 34
>>
>> Hence, indeed, some rings were garbage collected.
>
> I am not sure this fixes the problem mentioned above, but I am still
> testing. It might reduce the memory used already, but I won't know for
> a while.
>
>> Cheers,
>> Simon
>
> As is this patch breaks badly:

Yep, see

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

- Robert

[mabshoff]

On May 3, 1:20 am, Robert Bradshaw <rober...@math.washington.edu>
wrote:

> On May 3, 2009, at 1:04 AM, mabshoff wrote:

<SNIP>

> > As is this patch breaks badly:
>
> Yep, see
>
> http://trac.sagemath.org/sage_trac/ticket/3706
>
> - Robert

Do you mean we need an analog fix? I have always been mistrustful of
weak references, especially in Cython, since they seem to involve hard
to debug Heisenbugs in the code :)

Cheers,

Michael

[mabshoff]

FYI: Simon's patch doesn't make any difference for the amount of
memory used for


sage: while p<2^20:
....: p=next_prime(p)
....: g=FindGroupOrder(p,11)
....:

In both cases about 2046MB more were used after the loop :(. Simon:
Run

init_mem=get_memory_usage()

before you do anything, run the computation and then once you are done
run

get_memory_usage()-init_mem

and compare the different to before the patch. I assume Alex's script
he attached to the ticket does something similar.

It might be a good idea to create some constant at the very end of
all.py that sets the initial amount of memory used so it is trivially
invokable because once you have hit a leak it is too late anyway :(

Not that the patch will fix some other issues, but the picture seems
to be more complex. I am sure we will get to the bottom of it :).

Cheers,

Michael

[Alex Ghitza]

On Sun, May 3, 2009 at 6:33 PM, mabshoff
<Michael...@mathematik.uni-dortmund.de> wrote:
>
> FYI: Simon's patch doesn't make any difference for the amount of
> memory used for
>
>
> sage: while p<2^20:
> ....:     p=next_prime(p)
> ....:     g=FindGroupOrder(p,11)
> ....:
>
> In both cases about 2046MB more were used after the loop :(.

Note however that running just a loop where the polynomial ring is
created does benefit from the change in Simon's patch: doing it for p
up to 2^17, this used to eat up 272MB, and after the patch it's only
taking 10MB.

I did experience the same phenomenon that Michael is describing when
running loops where elliptic curves (or even just plane curves) are
created. The patch seems to have no effect on those loops.

[simon...@uni-jena.de]

Hi Michael,

On 3 Mai, 10:04, mabshoff <Michael.Absh...@mathematik.uni-dortmund.de>
wrote:

> Is there any reason you opened a new ticket and did not use #5949 as
> mentioned above?

Yes. I did the patch before breakfast :)

If it is true that there is a problem with weak references in Cython
(what are they?) then there will be no chance for my patch to work, I
guess.

Best regards,
Simon

[mabshoff]

On May 3, 2:00 am, simon.k...@uni-jena.de wrote:
> Hi Michael,
>
> On 3 Mai, 10:04, mabshoff <Michael.Absh...@mathematik.uni-dortmund.de>
> wrote:
>
> > Is there any reason you opened a new ticket and did not use #5949 as
> > mentioned above?
>
> Yes. I did the patch before breakfast :)

:)

Overall it seems to be good luck that this is on another ticket since
it seems to be part of the problem, but doesn't fix the issue in
question at the start of the thread.

> If it is true that there is a problem with weak references in Cython
> (what are they?) then there will be no chance for my patch to work, I
> guess.

It can be fixed - see the ticket RobertWB mentioned earlier in this
thread.

> Best regards,
>    Simon

Cheers,

Michael

[mabshoff]

On May 3, 1:54 am, Alex Ghitza <aghi...@gmail.com> wrote:
> On Sun, May 3, 2009 at 6:33 PM, mabshoff
>

> <Michael.Absh...@mathematik.uni-dortmund.de> wrote:
>
> > FYI: Simon's patch doesn't make any difference for the amount of
> > memory used for
>
> > sage: while p<2^20:
> > ....:     p=next_prime(p)
> > ....:     g=FindGroupOrder(p,11)
> > ....:
>
> > In both cases about 2046MB more were used after the loop :(.
>
> Note however that running just a loop where the polynomial ring is
> created does benefit from the change in Simon's patch: doing it for p
> up to 2^17, this used to eat up 272MB, and after the patch it's only
> taking 10MB.

Cool, that matches with my experience, so it is part of the solution
to get the initial problem fixed.

> I did experience the same phenomenon that Michael is describing when
> running loops where elliptic curves (or even just plane curves) are
> created.  The patch seems to have no effect on those loops.

Yep. I would not be surprised if Simon's patch just exposes the next
problem - oh well, fixing leaks is an eternal game of whack-a-mole ;).
I am all in favor on adding some looping tests that measure memory
consumption and that would fail if repeated runs start to leak memory
(again) once they are fixed.

> Alex

Cheers,

Michael