How to deal with GAP's machine dependent random generator?
Simon King
I almost accomplished to port my group cohomology spkg to T2.
As a side note: It was necessary to rename two functions in the
underlying C code, since on T2 (but nowhere else) apparently these
functions got confused with functions of the same name from pari. This
also turned out to be a problem for #1396.
The remaining problem concerns the random generator of GAP. It behaves
differently on T2 than on any other machine that I have access to.
On my desktop computer, on sage.math and on bsd.math (the latter both
in 64 and 32 bit mode), I get:
gap> Reset(GlobalMersenneTwister);; List([1..10],i->Random(1,100000));
[ 85758, 234, 18408, 12676, 2507, 38166, 24692, 15013, 4755, 94043 ]
On T2, I get:
gap> Reset(GlobalMersenneTwister);; List([1..10],i->Random(1,100000));
[ 96343, 43659, 75071, 1314, 40740, 69892, 10586, 31272, 96336,
26794 ]
It seems that due to this problem, my package may compute a different
(though isomorphic) ring presentation for a cohomology ring on
different machines.
But in some tests, I create elements with certain properties, and for
that purpose I need a computationally unique ring presentation, not
only for testing the result but for setting up the example. So, even
marking some test as random would not help.
What do you think is the best way to cope with that problem? Is there
a painless way to duplicate some doctests and tell the test script
that one version is for Solaris and the other is for everything else?
Best regards,
Simon
Dr. David Kirkby
> Hi!
>
> I almost accomplished to port my group cohomology spkg to T2.
>
> As a side note: It was necessary to rename two functions in the
> underlying C code, since on T2 (but nowhere else) apparently these
> functions got confused with functions of the same name from pari. This
> also turned out to be a problem for #1396.
>
> The remaining problem concerns the random generator of GAP. It behaves
> differently on T2 than on any other machine that I have access to.
>
> On my desktop computer, on sage.math and on bsd.math (the latter both
> in 64 and 32 bit mode), I get:
> gap> Reset(GlobalMersenneTwister);; List([1..10],i->Random(1,100000));
> [ 85758, 234, 18408, 12676, 2507, 38166, 24692, 15013, 4755, 94043 ]
>
> On T2, I get:
> gap> Reset(GlobalMersenneTwister);; List([1..10],i->Random(1,100000));
> [ 96343, 43659, 75071, 1314, 40740, 69892, 10586, 31272, 96336,
> 26794 ]
Do you need to use Random numbers? Can't you just pick a set of integers -
either those from t2, sage.math, or the lottery numbers, and use them instead?
Dave
Simon King
On Sep 6, 1:02 pm, "Dr. David Kirkby" <david.kir...@onetel.net> wrote:
> ...
> either those from t2, sage.math, or the lottery numbers, and use them instead?
way to show that there is a problem.
The randomness is only used internally. Namely, in my cohomology
programs happens the following:
(1) Input is some group.
(2) Internally, some additional data are computed, using GAP
functions such as IsomorphismGroups, SylowSubgroup etc. Some of these
functions use a randomized algorithm, and according to gap-support,
GAP does not provide a deterministic version of these algorithms.
(3) The ring presentation will depend on the internal data computed
in (2).
(4) In order to write tests (not only to evaluate test results), I
need a computationally unique machine independent ring presentation.
As the example shows, in order to achieve (4) it does not suffice to
reset GAP's random generator, since it is machine dependent.
Best regards,
Simon
Simon King
On Sep 6, 1:21 pm, Simon King <simon.k...@nuigalway.ie> wrote:
> > ...
> > Do you need to use Random numbers? Can't you just pick a set of integers -
> > either those from t2, sage.math, or the lottery numbers, and use them instead?
>
> I don't need these specific numbers. The example was just a concise
> way to show that there is a problem.
gap> S := SymmetricGroup(6);;
gap> G := SylowSubgroup(S,2);
Group([ (1,2), (3,4), (1,3)(2,4), (5,6) ])
gap> Id := IdGroup(G);;
and now, ``IsomorphismGroups(SmallGroup(Id),G);`` has the same result
on sage.math and bsd.math and my desktop computer, but a different
result on t2.
Cheers,
Simon
Simon King
> a painless way to duplicate some doctests and tell the test script
> that one version is for Solaris and the other is for everything else?
I.e.:
The spkg is distributed so that the tests work on sage.math, bsd.math,
hopefully PPC etc.
In addition, it contains a doctest patch: spkg-install verifies
whether we are on a Solaris OS, and if it is the case, the patch is
applied so that the different-but-isomorphic ring structures for
Solaris appear in the tests. With the patch, tests pass on t2.math as
well.
Do people think this is an acceptable solution?
Cheers,
Simon
Mike Hansen
> The remaining problem concerns the random generator of GAP. It behaves
> differently on T2 than on any other machine that I have access to.
My guess is that it is due to the endianness of the machine. For
example, check out this snippet from sage.misc.randstate:
# GAP's random number generator initialization
# (in integer.c, in FuncInitRandomMT) takes its
# seed as a string, then converts this string into
# an array of 32-bit integers just by casting the
# pointer. Thus, the result depends on the
# endianness of the machine. As a workaround, we
# swap the bytes in the string ourselves, so that
# GAP always gets the same array of integers.
If you're using GAP from withing Sage, you might just want to do
current_randstate().set_seed_gap() to initialize the random number
generators. Note that both of the following are reset:
prev_mersenne_seed = gap.Reset(gap.GlobalMersenneTwister,
mersenne_seed)
prev_classic_seed = gap.Reset(gap.GlobalRandomSource, classic_seed)
--Mike
Nils Bruin
> It seems that due to this problem, my package may compute a different
> (though isomorphic) ring presentation for a cohomology ring on
> different machines.
>
> But in some tests, I create elements with certain properties, and for
> that purpose I need a computationally unique ring presentation, not
> only for testing the result but for setting up the example. So, even
> marking some test as random would not help.
deterministic algorithm, they are going to be very fragile. They don't
just test the correctness of your code but some extra arbitrary
behaviour as well. If in the future, GAP would change one of its
internal algorithms you depend on, your test would almost certainly
break.
One possible way to handle this is by rewriting your tests in such a
way that you test the properties that matter (and should not depend on
some underlying implementation), in much the same way how set equality
is better tested by asking sage "are these sets equal?" than comparing
its string representation with a precomputed one.
In your case it sounds like you'd have to delve a bit deeper, because
you are currently making assumptions about the representation of the
ring to even get your example started. However, when you constructed
the example, you probably found your element with certain properties
by doing some computations and tests in the representation that sage
gave you. You could make those computations part of your test.
I realize that especially in cohomology rings, there is so much choice
in representation that making all those computations part of the test
might be unreasonable for a test-writer (I suspect almost certainly
not for testing time), so it's a judgement call for you: "Is a truly
robust test worth the effort from the test writer". If you decide
"no", you may want to mark the test with the ambiguity in the test
together with the reason why you left the ambiguity in.
Keep in mind that where subtle changes in behaviour of software occur
between platforms, there may well be similar changes in behaviour of
that software between different versions on the same platform.
Resetting random seeds is useful for reproducing errors. Less so for
reproducing test results.
Dima Pasechnik
endianness.
(I can reproduce this behaviour on a MacOSX PPC machine)
It can also be noted that apparently GAP uses an old Mersenne Twister
code, whereas it
is known to be buggy in exactly this way (endianness!), and fixes
exist...
If there is a pseudorandom numbers guru that can fix the GAP code in
GAPROOT/src/integer.c quickly, it would be great,
but it could be a while otherwise...
Using an external pseudorandom generator does not seem to be feasible,
as it's not going to be used
by GAP internals without much hacking.
Dima
On Sep 7, 1:09 am, Mike Hansen <mhan...@gmail.com> wrote:
Simon King
On 6 Sep., 19:44, Nils Bruin <nbr...@sfu.ca> wrote:
> If your tests depend on the choice of random numbers in a non-
> deterministic algorithm, they are going to be very fragile.
isomorphism between two finite groups" etc. is essential to the
computations. This is currently only available (vie GAP) by randomized
algorithms. So, I don't think I can easily change it, and have to
somehow cope with the instabilities (what I currently try to do).
> One possible way to handle this is by rewriting your tests in such a
> way that you test the properties that matter (and should not depend on
... yes, such as testing against the Hilbert-Poincaré series and other
invariants.
Sure. Probably I should use those kind of tests more often.
> However, when you constructed
> the example, you probably found your element with certain properties
> by doing some computations and tests in the representation that sage
> gave you. You could make those computations part of your test.
against theorems.
E.g., the program computes a homogeneous system of parameters, and
the test then demonstrates (as an example for the user) how to
conclude from properties of the parameters that the ring presentation
is complete. In these cases I often also show how the parameters look
like, but it would not harm to mark this part of the test as "random".
More difficult are the "easy" tests: Testing __mul__, __cmp__ etc,
because here one has no theorems to test against.
> Resetting random seeds is useful for reproducing errors. Less so for
> reproducing test results.
Cheers,
Simon
Simon King
On 6 Sep., 19:47, Dima Pasechnik <dimp...@gmail.com> wrote:
> On the GAP support list we reached the conclusion that it is due to
> endianness.
This is actually a little ironic: At Sage Days XX in San Diego, a
couple of years ago, both Martin Albrecht and Michael Abshoff said
that probably my wrapper for C-MeatAxe (which is part of the package)
would not work in different endianness, due to my seemingly flawed way
of pickling matrix data.
Now, it turns out that this is no problem. I can read MeatAxe matrices
on t2.math that I saved on sage.math.
On sage.math:
sage: from pGroupCohomology.mtx import MTX
sage: M = MTX(27,[[randint(0,26) for _ in range(5)] for _ in
range(5)])
sage: print M
[17 9 14 0 10]
[15 14 19 3 19]
[16 3 19 23 8]
[ 7 24 3 2 24]
[19 3 20 1 2]
sage: save(M, 'matrixtest')
On t2:
sage: M = load('matrixtest.sobj')
sage: print M
[17 9 14 0 10]
[15 14 19 3 19]
[16 3 19 23 8]
[ 7 24 3 2 24]
[19 3 20 1 2]
But endianness strikes back in an unexpected way...
Cheers,
Simon
Dr. David Kirkby
> Hi Dima!
>
> On 6 Sep., 19:47, Dima Pasechnik<dimp...@gmail.com> wrote:
>> On the GAP support list we reached the conclusion that it is due to
>> endianness.
>
> Is the endianness really different on sage.math and t2.math?
Yes. x86 are little endian. SPARC is big endian.
I can try the code if you want on my OpenSolaris system. That uses an Intel Xeon
processor, so is little endian.
I used to work on a 24-bit machine. Either the most (or least), my memory fails
me, was in the middle!
I would be very weary of any random number generator that claims to be a good
source of random numbers if the output differs by platform or compilation mode.
Someone would need to check the numbers for randomless on these different systems.
There are a number of statistical tests for pseudo random number generators -
but I expect you know that anyway.
> Cheers,
> Simon
>
Dave
Simon King
On 6 Sep., 22:29, "Dr. David Kirkby" <david.kir...@onetel.net> wrote:
> > Is the endianness really different on sage.math and t2.math?
>
> Yes. x86 are little endian. SPARC is big endian.
> I can try the code if you want on my OpenSolaris system. That uses an Intel Xeon
> processor, so is little endian.
> > But endianness strikes back in an unexpected way...
>
> I used to work on a 24-bit machine. Either the most (or least), my memory fails
> me, was in the middle!
Best regards,,
Simon
Dr. David Kirkby
> Hi Dave!
>
> On 6 Sep., 22:29, "Dr. David Kirkby"<david.kir...@onetel.net> wrote:
>>> Is the endianness really different on sage.math and t2.math?
>>
>> Yes. x86 are little endian. SPARC is big endian.
>
> Cool! Now I am really proud that my code mostly works on t2.
>
>> I can try the code if you want on my OpenSolaris system. That uses an Intel Xeon
>> processor, so is little endian.
>
> Reset(GlobalMersenneTwister);; List([1..10],i->Random(1,100000));
This is on a Sun Ultra 27 running OpenSolaris, with a Xeon box.
sage -t -long devel/sage/sage/quadratic_forms/quadratic_form__siegel_product.py
[164.1 s]
----------------------------------------------------------------------
All tests passed!
Total time for all tests: 1825.6 seconds
drkirkby@hawk:~/noatlas/sage-4.5.3.rc0$ ./sage -gap
######### ###### ########### ###
############# ###### ############ ####
############## ######## ############# #####
############### ######## ##### ###### #####
###### # ######### ##### ##### ######
###### ########## ##### ##### #######
##### ##### #### ##### ###### ########
#### ##### ##### ############# ### ####
##### ####### #### #### ########### #### ####
##### ####### ##### ##### ###### #### ####
##### ####### ##### ##### ##### #############
##### ##### ################ ##### #############
###### ##### ################ ##### #############
################ ################## ##### ####
############### ##### ##### ##### ####
############# ##### ##### ##### ####
######### ##### ##### ##### ####
Information at: http://www.gap-system.org
Try '?help' for help. See also '?copyright' and '?authors'
Loading the library. Please be patient, this may take a while.
GAP4, Version: 4.4.12 of 17-Dec-2008, i386-pc-solaris2.11-gcc
gap> Reset(GlobalMersenneTwister);; List([1..10],i->Random(1,100000));
[ 85758, 234, 18408, 12676, 2507, 38166, 24692, 15013, 4755, 94043 ]
gap>
Simon King
On 6 Sep., 23:53, "Dr. David Kirkby" <david.kir...@onetel.net> wrote:
> ...
>
> This is on a Sun Ultra 27 running OpenSolaris, with a Xeon box.
>
> [ 85758, 234, 18408, 12676, 2507, 38166, 24692, 15013, 4755, 94043 ]
So indeed, it is not "Solaris vs. non-Solaris" but "little vs. big
endian", as the GAP people suspected.
Cheers,
Simon
Mike Hansen
> So indeed, it is not "Solaris vs. non-Solaris" but "little vs. big
> endian", as the GAP people suspected.
This is what I said in my first response and can be taken care of with
current_randstate().set_seed_gap().
Core2:
sage: set_random_seed(100)
sage: current_randstate().set_seed_gap()
sage: gap.eval('List([1..10],i->Random(1,100000));')
'[ 80761, 80557, 8462, 25730, 32000, 92371, 85402, 52347, 17181, 2889 ]'
t2:
sage: set_random_seed(100)
sage: current_randstate().set_seed_gap()
sage: gap.eval('List([1..10],i->Random(1,100000));')
'[ 80761, 80557, 8462, 25730, 32000, 92371, 85402, 52347, 17181, 2889 ]'
--Mike
Dima Pasechnik
sage: set_random_seed(100)
sage: current_randstate().set_seed_gap()
sage: gap.eval('List([1..10],i->Random(1,100000));')
'[ 80761, 80557, 8462, 25730, 32000, 92371, 85402, 52347, 17181,
2889 ]'
GAP4, Version: 4.4.12 of 17-Dec-2008, powerpc-apple-darwin9.8.0-gcc
gap> List([1..10],i->Random(1,100000));
26794 ]
I wonder how does GAP pick up the seed? (Don't even know where to look
in the code...)
On the other hand, how would one call
current_randstate().set_seed_gap() from GAP?
(that's the thing that would be needed to fix GAP's misbehaviour...)
Dima
On Sep 7, 7:05 am, Mike Hansen <mhan...@gmail.com> wrote:
Mike Hansen
> I wonder how does GAP pick up the seed? (Don't even know where to look
> in the code...)
It is in sage/misc/randstate.pyx.
> On the other hand, how would one call
> current_randstate().set_seed_gap() from GAP?
You can't call it from GAP, but it just does the following calls
prev_mersenne_seed = gap.Reset(gap.GlobalMersenneTwister,
mersenne_seed)
prev_classic_seed = gap.Reset(gap.GlobalRandomSource, classic_seed)
which are just GAP calls. If the machine is little endian, then
classic_seed and mersenne_seed are the same. If the machine is big
endian, then the following following comment applies:
# GAP's random number generator initialization
# (in integer.c, in FuncInitRandomMT) takes its
# seed as a string, then converts this string into
# an array of 32-bit integers just by casting the
# pointer. Thus, the result depends on the
# endianness of the machine. As a workaround, we
# swap the bytes in the string ourselves, so that
# GAP always gets the same array of integers.
Thus, the following is done
seed = str(seed)
new_seed = ''
while len(seed) >= 4:
new_seed += seed[3::-1]
seed = seed[4:]
seed = '"' + new_seed + '"'
mersenne_seed = seed
--Mike
Simon King
On 7 Sep., 01:05, Mike Hansen <mhan...@gmail.com> wrote:
I did not expect that it is so carefully implemented - I thought it is
just calling the GAP commands, without a special case for endianness.
So, I will try this. Probably I have to rewrite the tests anyway (the
above sequence is different from the one that I get with
Reset(GlobalMersenneTwister), but at least it probably means that
eventually the ring structures will be computationally unique on all
platforms.
Best regards,
Simon
Simon King
> That's VERY cool!
... modulo one critical comment: In principal (although it should
rarely occur in "real life") one can have several instances of an
interface. So, I think that current_randstate().set_seed_gap() should
accept an optional argument (namely an instance of GAP), and only if
it is not provided then the default sage.interfaces.gap should be
chosen.
Do people feel I should modify it accordingly?
Cheers,
Simon
Dima Pasechnik
I've made a small patch to GAP that fixes the problem.
(inspired by the corresponding Sage Python code :))
With this patch I get on Sparc:
sage: gap_console()
GAP4, Version: 4.4.12 of 17-Dec-2008, sparc-sun-solaris2.10-gcc
[ 85758, 234, 18408, 12676, 2507, 38166, 24692, 15013, 4755, 94043 ]
the byte order in the appropriate integers:
The following is the diff for GAPROOT/src/integer.c
(the context diff is attached)
gap4r4/src$ diff integer.c.old integer.c
3120c3120
< UInt4 *mt, *init_key, key_length, i, j, k, N=624;
---
> UInt4 *mt, *init_key, key_length, i, j, k, N=624, temp;
3140a3141,3146
> #ifdef WORDS_BIGENDIAN /* to make the behaviour uniform across platforms */
> temp = ((init_key[j] & 0xFF) << 24) | ((init_key[j] & 0xFF00) << 8)
> | ((init_key[j] >> 8) & 0xFF00) | ((init_key[j] >> 24) & 0xFF);
> #else
> temp = init_key[j];
> #endif
3142c3148
< + init_key[j] + j;
---
> + temp + j;
It would be nice to put this in the upcoming GAP release...
====== the following is Sage-specific: ===
I have put the updated Sage package file gap-4.4.12.p4.spkg
(I didn't change the number nor commited changes yet) here:
http://boxen.math.washington.edu/home/dima/packages/gap-4.4.12.p4.spkg
(so you can try and see if this works for you)
To update sage with it correctly,
one should also remove the analogous fix in
SAGEROOT/devel/sage/sage/misc/randstate.pyx
(lines 687-703, starting from "if sys.byteorder == 'big':"
in set_seed_gap)
If you like, we can quickly make another update of gap spkg with this
fix (in sync with the fix for randstate.pyx)
====================
Best,
Dima
Dima Pasechnik
http://boxen.math.washington.edu/home/dima/packages/gap-4.4.12.p5.spkg
with the proper patches, bells, etc., and
will make a ticket to update the things shortly.
Dima Pasechnik
Please review!
On Sep 7, 8:56 pm, Dima Pasechnik <dimp...@gmail.com> wrote:
> I've madehttp://boxen.math.washington.edu/home/dima/packages/gap-4.4.12.p5.spkg
Peter Jeremy
>I would be very weary of any random number generator that claims to be a good
>source of random numbers if the output differs by platform or compilation mode.
It depends what you mean by "differs". For "real" random numbers, you
can only discuss different statistical properties of the sequences
because the actual sequencee of numbers cannot be predicted or repeated.
For pseudo-random sequences, you can additionally talk about the
actual sequence of numbers returned by the rng. The documentation for
most rng's talks about sequences being repeatable by calling an
initialisation function with the same seed but generally does not
state the scope of the repeatability - so actual sequences generated
by a specific seed could potentially vary on different platforms.
I would suggest that the bigger problem is that the GAP doctest relies
on a specific sequence of random numbers being used by GAP but does
not control or test the actual sequence of random numbers.
--
Peter Jeremy
Dima Pasechnik
On Sep 9, 3:12 am, Peter Jeremy <peterjer...@acm.org> wrote:
GAP, without the patch I proposed in the spkg in this thread (trac
#9867), uses a random seed in an endianness-dependent way.
The said patch fixes this.
Nobody seems to be willing to review the ticket, however, as if the
problem got resolved in a way that escapes me.
Best,
Dima
> --
> Peter Jeremy
>
> application_pgp-signature_part
> < 1KViewDownload
Simon King
On Sep 9, 2:44 pm, Dima Pasechnik <dimp...@gmail.com> wrote:
>...
> problem got resolved in a way that escapes me.
current_randstate().set_seed_gap(). But of course it should better be
solved directly in GAP (resp. the GAP spkg).
However, at the moment, I spend all my energy in (1) finalising the
next version of my cohomology package, and (2) the relocation of my
family (without the help of my wife, who is very sick). After (1) I'll
probably need a break from development. So, I'd appreciate if someone
else could review the ticket.
Best regards,
Simon
Dima Pasechnik
On Sep 9, 10:23 pm, Simon King <simon.k...@nuigalway.ie> wrote:
> Hi Dima!
>
> On Sep 9, 2:44 pm, Dima Pasechnik <dimp...@gmail.com> wrote:
>
> >...
> > Nobody seems to be willing to review the ticket, however, as if the
> > problem got resolved in a way that escapes me.
>
> The problem for my package did indeed resolve, due to
> current_randstate().set_seed_gap(). But of course it should better be
> solved directly in GAP (resp. the GAP spkg).
seed, so that
subsequent GAP commands make use of the correctly set seed?
>
> However, at the moment, I spend all my energy in (1) finalising the
> next version of my cohomology package, and (2) the relocation of my
> family (without the help of my wife, who is very sick). After (1) I'll
> probably need a break from development. So, I'd appreciate if someone
> else could review the ticket.
Dima
>
> Best regards,
> Simon
Mike Hansen
> Does current_randstate().set_seed_gap() actually sets GAP's random
> seed, so that
> subsequent GAP commands make use of the correctly set seed?
Yep.
--Mike
Dima Pasechnik
On Sep 11, 2:36 pm, Mike Hansen <mhan...@gmail.com> wrote:
> > Does current_randstate().set_seed_gap() actually sets GAP's random
> > seed, so that
> > subsequent GAP commands make use of the correctly set seed?
>
> Yep.
ticket duly reflects),
so this will have to be dealt with when the new GAP is released
anyway.
Dima
>
> --Mike