Comparison of Integer & numpy.int32 fails... but not always?

[Maxim]

Hello everyone,

I'd like to know if it is required practice to cast the types of the
objects being compared or if Sage should be able to deal with it
automatically?

Because I have a problem where I need to explicitely cast an
numpy.int32 vector unit for the comparison with a Sage Integer to be
reliable. It would sometimes see a '1' Integer not being the same as a
'1' numpy.int32 value.

This behaviour _cannot_ be reproduced with a minimal example such as:

import numpy

if numpy.int32(1) != sage.rings.integer.Integer(1):
print("Should not go here...")


If someone wants to dig to the root of the problem, I can gladly
provide my complete program.

[Maxim]

Well, I'd like to dig the root of the problem ;). So I'm listing the
pertinent code.

Also if you see performance improvement possible optimization to my
code, they are very welcome, as this
thing is rather slow.

Abstract: Requirement was to build a pseudo random bit generator using
linear congruence,
and also generate gaussian distributed "noise" using Teichroew's
method to simulate transmission
of non modulated binary data in a noisy channel. Note: all error
control code was removed to simplify
the listing (and my comments translation job ;)). It isn't really
needed anyway, unless you play with the parameters of
my functions.

## START OF CODE ##
# We use numpy vectors because of ticket #10262.
import numpy
"""
Random integers generator using linear congruence
Prototype: LCG(n,x=1,a=16807,c=0,m=2^31-1)
n: number of pseudo random integers we want to generate
x: initial seed
a: multiplier (a < m)
c: increment( c < m)
m: module
"""
def LCG(n,x=1,a=16807,c=0,m=2^31-1):
rand_int_list=[]
for i in range(n):
x = (a*x+c) % m
rand_int_list.append(x)
return numpy.array(rand_int_list)

# Normalized ([0,1[) pseudo random integers
def LCG_norm(n,x=1,a=16807,c=0,m=2^31-1):
return LCG(n,x,a,c,)/float(m-1)

# Pseudo random bit generator (we keep a single bit out of a LCG
generated integer)
# p: weight of the bit we want to isolate
def LCG_bits(n,p=3,x=1,a=16807,c=0,m=2^31-1):
int_vect = LCG(n,x,a,c,m)
for i in range(len(int_vect)):
int_vect[i] = (int_vect[i] & 2^p)/2^p
return int_vect

# Gaussian distributed integer numbers generator using Teichroew's
method
def rand_gauss(n,moy=0,sigma=1):
LCG_uni_vector=LCG_norm(n*12)
n_k=[]
for k in range(n):
somme=0
for i in range(k*12,k*12+12):
somme+=LCG_uni_vector[i]
R=(somme-6)/4
a1=3.949846138; a3=0.252408784
a5=0.076542912; a7=0.008355968
a9=0.029899776
X_T = ((((a9*R^2 + a7)*R^2 + a5)*R^2 + a3)*R^2 + a1)*R
n_k.append(X_T*sigma + moy)
return numpy.array(n_k)

# Transmission probability of error simulation using Monte Carlo's
method.
# The detector is a simple decision rule based on the level of the
signal + noise.
def erreur_transmission(n,A=1,sigma=1):
vecteur_bits=LCG_bits(n)
vecteur_bruit=rand_gauss(n,0,sigma)
seuil=A/2
nb_erreurs=0
for k in range(n):
y_k = A*vecteur_bits[k] + vecteur_bruit[k]
r_k = 1 if (y_k > seuil) else 0
# EXPLICIT CASTING WAS FOUND TO BE REQUIRED HERE... ELSE
INCONSISTENT RESULTS
if Integer(r_k) != Integer(vecteur_bits[k]):
nb_erreurs = nb_erreurs+1
P_erreur = nb_erreurs/n
return float(P_erreur)
## END OF CODE ##

Now you can try the program by issuing:
$> erreur_transmission(100)

You should get a value of: 0.25999999999999995 (this correspond to the
probability of error calculated with a length of 100 bits.
sigma & A values were left to there default, which gives a SNR of
~-6dB.)

Now the interesting experiment: in the last function, under the
capitalized comment, line, remove the Integer casting like this:
if r_k != vecteur_bits[k]:

Rerun the "erreur_transmission(100)" command.

You should get a value of: 0.59999999999999998.

And this is what I see as an inconsistent and buggy behaviour. Or is
there a reason for this?

[Maxim]

Google groups long lines wrapping breaks my code. Check this text
version if you don't want to bother refactoring it:

http://dl.dropbox.com/u/176092/comparison_test.sage

[D. S. McNeil]

On Wed, Feb 16, 2011 at 5:53 AM, Maxim <maxim.c...@gmail.com> wrote:

> It would sometimes see a '1' Integer not being the same as a
> '1' numpy.int32 value.
>
> This behaviour _cannot_ be reproduced with a minimal example such as:
>
> import numpy
>
> if numpy.int32(1) != sage.rings.integer.Integer(1):
>    print("Should not go here...")

I think it can, but only in the other order:

sage: if Integer(1) != numpy.int32(1):
....: print("should not go here")
....:
should not go here

I don't think it's a vector/non-vector issue. Probably following
whatever cmp approach Integers use will eventually bring
enlightenment..

[As for speed, you seem to be using a lot of 'for i in range(n)'
loops, some of which can be pushed into vector operations.]

Doug

--
Department of Earth Sciences
University of Hong Kong