numpy-discussion — Discussion list for all users of Numerical Python

[Numpy-discussion] Re: numarray: RandomArray2

[Jochen <jo...@jo...>]

On Sun, 18 Aug 2002 08:13:18 -0400 Todd Miller wrote:

Todd> Jochen Küpper wrote:

>> Looking at RandomArray2 I realized that the functions don't do any
>> argument checking. Shouldn't uniform for example check at least
>> whether ,----
>> | minimum != maximum
>> `----
>> or even make sure that maximum > minimum?  (Although the result is
>> probably "ok" if it isn't.)  Something like
>> ,----
>> | def uniform(minimum, maximum, shape=[]):
>> |     """Return array of random Floats in interval ]minimum, maximum[
>> |     with shape |shape|.
>> |     """
>> |     if maximum <= minimum:
>> |         raise ValueError
>> |     return minimum + (maximum-minimum)*random(shape)
>> `----
>> 
Todd> I am +0 on this. The parameter order emulates "range". 

Not exactly, currently:
,----
| >>> range(2, -4)
| []
| >>> ran.uniform(2, -4)
| -2.2346744537353516
`----
That is, if max < min range returns an empty list, whereas uniform
comes up with the same result as for (-4, 2) (well, "mirrored").

Also note the "strange" docstring of range:
,----
| range(...)
|     range([start,] stop[, step]) -> list of integers
`----
pointing straight toward its behavior.

Todd> While it does look like it will compute the wrong answer if
Todd> called incorrectly, that seems unlikely to me.

Well, if max<min it "mirrors" the result inside the range. This is no
problem as long as it is done "consistently". I don't have enough
knowledge of RNG's to understand whether it is a problem (with respect
to randomness) when calls with the right and wrong ordering of min and
max are mixed.

Thinking about the case min==max I must say it's a very wasted
function call, but no actually big deal:
,----
| >>> import RandomArray2 as ran
| >>> ran.uniform(1, 1)
| 1.0
`----

So well, maybe someone with insight into RNG's can comment on the,
mirroring issue?

>> Moreover there are some inconsistencies between functions, i.e.:
>> ,----
>> | def randint(minimum, maximum=None, shape=[]):
>> `----
>> ,----
>> | def random_integers(maximum, minimum=1, shape=[]):
>> `----
 
Todd> It appears to me that the parameter order of randint again
Todd> emulates "range". The fact that random_integers is not
Todd> consistent with randint seem OK to me because random_integers
Todd> appears to have been written expressly to tailor the calling
Todd> sequence of randint.

Hmm,

v.s. 

Yes, initially I wondered why there are two functions at all. This
explanation sounds like "we wanted some help in confusing the
users"  :))

Todd> Because randint re-defines its parameters depending on whether 1
Todd> or 2 range values are used, as does range, I don't think there
Todd> is a completely consistent way to do this. Either we're "wrong"
Todd> for the 1 parameter case or the 2 parameter case. The way it is
Todd> specified now seems simplest to me, with "minimum" preceding
Todd> "maximum", even though it is not strictly the correct name for
Todd> the 1 parameter case.

What's about 
,----
| uniform(limit1, limit2, shape=[])
`----
and then range-like behaviour?

But then, what do we return on 
,----
| uniform(2, -4)
`----
??? 

To make it compatible with range, it should be an empty list, no?

Your-even-more-confused-ly's,
Jochen
-- 
Einigkeit und Recht und Freiheit                http://www.Jochen-Kuepper.de
    Liberté, Égalité, Fraternité                GnuPG key: 44BCCD8E
        Sex, drugs and rock-n-roll

[Numpy-discussion] [ANNOUNCE] Pyfort 8.2 and video tutorials for Python, Pyfort

[Paul F D. <pa...@pf...>]

Pyfortran 8.2 is available. The setup bug for Windows is fixed and a -n
option has been added to just generate glue code and stop.

Video tutorials are now available at http://esg.llnl.gov/cdat/ (go to
documentation and tutorials).  There are instructions online about how
to set up your computer to play these videos. The Linux instructions may
or may not work for you; for the moment we recommend viewing from
Windows.

These tutorials are mostly for CDAT users, but there are two for Pyfort
and one short "Introduction to Python". Our most excellent summer
student from USC, Tasha Drew, got everything working.

 

Re: [Numpy-discussion] QR algorithms in Python

[Nils W. <nw...@me...>]

James G Analytis schrieb:
> 
> Dear Pearu,
> SciPy's linalg modules are just what I need. However, I get the following
> problem on importing SciPy:
> 
ATLAS is also required.

http://math-atlas.sourceforge.net/

Nils

> [analytis@toomey analytis]$ python2
> Python 2.2 (#1, Apr 12 2002, 15:29:57)
> [GCC 2.96 20000731 (Red Hat Linux 7.2 2.96-109)] on linux2
> Type "help", "copyright", "credits" or "license" for more information.
> >>> import scipy
> exceptions.ImportError:
> /usr/lib/python2.2/site-packages/scipy/linalg/flapack.so     : undefined
> symbol: sgesdd_
> exceptions.ImportError:
> /usr/lib/python2.2/site-packages/scipy/linalg/_flinalg.s     o: undefined
> symbol: dlaswp_
> Traceback (most recent call last):
>   File "<stdin>", line 1, in ?
>   File "/tmp/SciPyTest/linux2/lib/python2.2/site-packages/scipy/__init__.py",
> li     ne 42, in ?
>   File "/usr/lib/python2.2/site-packages/scipy/special/__init__.py", line 328,
> i     n ?
>     import orthogonal
>   File "/usr/lib/python2.2/site-packages/scipy/special/orthogonal.py", line
> 59,      in ?
>     from scipy.linalg import eig
>   File
> "/tmp/SciPyTest/linux2/lib/python2.2/site-packages/scipy/linalg/__init__.
> py", line 40, in ?
>   File
> "/tmp/SciPyTest/linux2/lib/python2.2/site-packages/scipy/linalg/basic.py"
> , line 17, in ?
> ImportError: /usr/lib/python2.2/site-packages/scipy/linalg/calc_lwork.so:
> undefi     ned symbol: ieeeck_
> >>>
> 
> I am running standard Red Hat 7.3 with a smp kernel, on a dual-processor PIII.
> Is there anything I should watch out for (compiler issues perhaps)?
> I appreciate the help!
> Cheers,
> James
> 
> On Monday 12 Aug 2002 9:04 am, Pearu Peterson wrote:
> > On Mon, 12 Aug 2002, James G Analytis wrote:
> > > Hi,
> > > I'm looking for a QR algorithm that can decompose complex matrices. The
> > > LinearAlgebra package from Numeric doesn't seem to have it. SciPy's
> > > linalg package is still on the way, so I can't use that either.
> >
> > Try SciPy's linalg package again (from CVS), the QR functions is working
> >
> > now:
> > >>> from scipy import linalg
> > >>> from Numeric import dot
> > >>> q,r=linalg.qr([[1,2],[3j,4]])
> > >>> print q
> >
> > [[-0.31622777+0.j         -0.78935222-0.52623481j]
> >  [ 0.        -0.9486833j  -0.1754116 +0.26311741j]]
> >
> > >>> print r
> >
> > [[-3.16227766+0.j         -0.63245553+3.79473319j]
> >  [ 0.        +0.j         -2.28035085+0.j        ]]
> >
> > >>> print dot(q,r)
> >
> > [[ 1. +0.00000000e+00j  2. -8.71264866e-16j]
> >  [ 0. +3.00000000e+00j  4. +3.09051829e-16j]]
> >
> > Regards,
> >       Pearu
> 
> -------------------------------------------------------
> This sf.net email is sponsored by: Dice - The leading online job board
> for high-tech professionals. Search and apply for tech jobs today!
> http://seeker.dice.com/seeker.epl?rel_code1
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] QR algorithms in Python

[James G A. <jam...@ph...>]

Dear Pearu,
SciPy's linalg modules are just what I need. However, I get the following=
=20
problem on importing SciPy:

[analytis@toomey analytis]$ python2
Python 2.2 (#1, Apr 12 2002, 15:29:57)
[GCC 2.96 20000731 (Red Hat Linux 7.2 2.96-109)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import scipy
exceptions.ImportError:=20
/usr/lib/python2.2/site-packages/scipy/linalg/flapack.so     : undefined=20
symbol: sgesdd_
exceptions.ImportError:=20
/usr/lib/python2.2/site-packages/scipy/linalg/_flinalg.s     o: undefined=
=20
symbol: dlaswp_
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File "/tmp/SciPyTest/linux2/lib/python2.2/site-packages/scipy/__init__.=
py",=20
li     ne 42, in ?
  File "/usr/lib/python2.2/site-packages/scipy/special/__init__.py", line=
 328,=20
i     n ?
    import orthogonal
  File "/usr/lib/python2.2/site-packages/scipy/special/orthogonal.py", li=
ne=20
59,      in ?
    from scipy.linalg import eig
  File=20
"/tmp/SciPyTest/linux2/lib/python2.2/site-packages/scipy/linalg/__init__.=
    =20
py", line 40, in ?
  File=20
"/tmp/SciPyTest/linux2/lib/python2.2/site-packages/scipy/linalg/basic.py"=
    =20
, line 17, in ?
ImportError: /usr/lib/python2.2/site-packages/scipy/linalg/calc_lwork.so:=
=20
undefi     ned symbol: ieeeck_
>>>

I am running standard Red Hat 7.3 with a smp kernel, on a dual-processor =
PIII.=20
Is there anything I should watch out for (compiler issues perhaps)?
I appreciate the help!
Cheers,
James

On Monday 12 Aug 2002 9:04 am, Pearu Peterson wrote:
> On Mon, 12 Aug 2002, James G Analytis wrote:
> > Hi,
> > I'm looking for a QR algorithm that can decompose complex matrices. T=
he
> > LinearAlgebra package from Numeric doesn't seem to have it. SciPy's
> > linalg package is still on the way, so I can't use that either.
>
> Try SciPy's linalg package again (from CVS), the QR functions is workin=
g
>
> now:
> >>> from scipy import linalg
> >>> from Numeric import dot
> >>> q,r=3Dlinalg.qr([[1,2],[3j,4]])
> >>> print q
>
> [[-0.31622777+0.j         -0.78935222-0.52623481j]
>  [ 0.        -0.9486833j  -0.1754116 +0.26311741j]]
>
> >>> print r
>
> [[-3.16227766+0.j         -0.63245553+3.79473319j]
>  [ 0.        +0.j         -2.28035085+0.j        ]]
>
> >>> print dot(q,r)
>
> [[ 1. +0.00000000e+00j  2. -8.71264866e-16j]
>  [ 0. +3.00000000e+00j  4. +3.09051829e-16j]]
>
> Regards,
> =09Pearu

Re: [Numpy-discussion] Scientific data - Reading files

[Chris B. <Chr...@no...>]

Nils Wagner wrote:
> Each file  (E1431A*.DAT) contains
> one column with real numbers. The contents of the i-th file should be
> stored in the i-th row of a data matrix A.

Are the files binary or ASCI ?

if ASCI, look in Scipy, there are a couple of modules for this kind of
thing, I doubt any of them are set up for multiple files, so you'll have
to loop through the file list.

If binary, you can do something like:

while row < NumRows:
	file = open(filename,'rb')
	A[row,:] = fromstring(file.read(),Float) # make sure you get the
typecode right!
	file.close()
	row += 1

WARNING: not the least bit tested.

-Chris



-- 
Christopher Barker, Ph.D.
Oceanographer
                                    		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chr...@no...

Re: [Numpy-discussion] QR algorithms in Python

[Konrad H. <hi...@cn...>]

James G Analytis <jam...@ph...> writes:

> I cannot compile PyLapack because I get this problem:
> 
> [analytis@toomey Modules]$ python2 compile.py
> Copy Misc/Makefile.pre.in from the Python distribution
> to this directory and try again.
> make: Makefile.pre.in: No such file or directory

The old module compilation procedure via Makefile.pre.in does not exist
any more in Python 2.2. If and when I find the time, I'll write a
distutils script for PyLapack.

Konrad.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: hi...@cn...
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------

Re: [Numpy-discussion] floating point exceptions

[Konrad H. <hi...@cn...>]

bryan cole <bry...@te...> writes:

> I'm trying to generate a gaussian distribution by:
> 
> Input = Numeric.arange(-1024,1024)/20.0
> Output = Numeric.exp(-(Input**2))
> 
> At the edges of the distribution the largish negative arguments cause
> the f.p. exception and prevent calculation of the entire array.
> 
> On windows the underflow errors are presumably rounded down to 0.0
> automatically. Is it possible to get this behaviour on linux/ix86?

Yes, rebuild the Python interpreter (not NumPy!) and add the compilation
option "-lieee". Alternatively, download and install the current
development version of Python from the CVS server (I haven't tried it,
but it is supposed to fix this problem as well).

Ultimately the problem is in the C standard library. On some platforms
it signals underflow as an overflow error, on others it doesn't. Python
thus has to get work around such problems in a platform-specific way.

Konrad.
-- 
-------------------------------------------------------------------------------
Konrad Hinsen                            | E-Mail: hi...@cn...
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------

Re: [Numpy-discussion] QR algorithms in Python

[Pearu P. <pe...@ce...>]

On Mon, 12 Aug 2002, James G Analytis wrote:

> Hi,
> I'm looking for a QR algorithm that can decompose complex matrices. The 
> LinearAlgebra package from Numeric doesn't seem to have it. SciPy's linalg 
> package is still on the way, so I can't use that either. 

Try SciPy's linalg package again (from CVS), the QR functions is working
now:

>>> from scipy import linalg
>>> from Numeric import dot
>>> q,r=linalg.qr([[1,2],[3j,4]])
>>> print q
[[-0.31622777+0.j         -0.78935222-0.52623481j]
 [ 0.        -0.9486833j  -0.1754116 +0.26311741j]]
>>> print r
[[-3.16227766+0.j         -0.63245553+3.79473319j]
 [ 0.        +0.j         -2.28035085+0.j        ]]
>>> print dot(q,r)
[[ 1. +0.00000000e+00j  2. -8.71264866e-16j]
 [ 0. +3.00000000e+00j  4. +3.09051829e-16j]]

Regards,
	Pearu

[Numpy-discussion] Scientific data - Reading files

[Nils W. <nw...@me...>]

Hi,

I have collected acceleration data by a multi-channel data acquisition
system. Each file  (E1431A*.DAT) contains
one column with real numbers. The contents of the i-th file should be
stored in the i-th row of a data matrix A.
How can I manage this task efficiently ? Are there suitable python
modules for this purpose ? 
A small example would be appreciated.

Nils

[Numpy-discussion] QR algorithms in Python

[James G A. <jam...@ph...>]

Hi,
I'm looking for a QR algorithm that can decompose complex matrices. The=20
LinearAlgebra package from Numeric doesn't seem to have it. SciPy's linal=
g=20
package is still on the way, so I can't use that either. Matfunc's QR=20
decomposition only works for real matrices. This problem came up two year=
s=20
ago on this same discussion panel, and Konrad Hinsen recommended that the=
=20
user download PyLapack from ftp://dirac.cnrs-orleans.fr/pub/
I cannot compile PyLapack because I get this problem:

[analytis@toomey Modules]$ python2 compile.py
Copy Misc/Makefile.pre.in from the Python distribution
to this directory and try again.
make: Makefile.pre.in: No such file or directory
make: *** No rule to make target `Makefile.pre.in'.  Stop.
make: *** No targets specified and no makefile found.  Stop.
[analytis@toomey Modules]$

I think the Makefile searched for is one from an earlier version of Numer=
ic=20
(the tar file is 2yrs old). Any suggestions?
Cheers,
James

Re: [Numpy-discussion] floating point exceptions

[bryan c. <bry...@te...>]

OK, I've worked out how to do this with MA.

I mask the Numeric array input to MA.exp() using MA.masked_outside( ...
,-700, 700) then convert the result back to a Numeric array using
MA.filled( ..., 0.0) then continue as normal.

I'm still interested if there's a more elegant approach.

Bryan

On Fri, 2002-08-09 at 16:52, bryan cole wrote:
> I've just ported a Numeric+Python program from windows to linux and I'm
> generating floating point exceptions in the linux version which did not
> occur in windows. 
> 
> I'm trying to generate a gaussian distribution by:
> 
> Input = Numeric.arange(-1024,1024)/20.0
> Output = Numeric.exp(-(Input**2))
> 
> At the edges of the distribution the largish negative arguments cause
> the f.p. exception and prevent calculation of the entire array.
> 
> On windows the underflow errors are presumably rounded down to 0.0
> automatically. Is it possible to get this behaviour on linux/ix86?
> 
> Anyone got any suggestions on a workarround for this? I see the numpy
> docs suggest using Masked Arrays, but since I don't know in advance
> which array items will generate the exceptions, I don't know which items
> to mask. Can the MA module generate an array mask based on f.p.
> exceptions from a calculation?
> 
> any suggestions would be much appreciated.
> 
> Bryan
> 
> 
> -- 
> Bryan Cole
> Teraview Ltd., 302-304 Cambridge Science Park, Milton Road, Cambridge
> CB4 0WG, United Kingdom.
> tel: +44 (1223) 435380 / 435386 (direct-dial) fax: +44 (1223) 435382
> 
> 
> 
> -------------------------------------------------------
> This sf.net email is sponsored by:ThinkGeek
> Welcome to geek heaven.
> http://thinkgeek.com/sf
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
-- 
Bryan Cole
Teraview Ltd., 302-304 Cambridge Science Park, Milton Road, Cambridge
CB4 0WG, United Kingdom.
tel: +44 (1223) 435380 / 435386 (direct-dial) fax: +44 (1223) 435382

[Numpy-discussion] floating point exceptions

[bryan c. <bry...@te...>]

I've just ported a Numeric+Python program from windows to linux and I'm
generating floating point exceptions in the linux version which did not
occur in windows. 

I'm trying to generate a gaussian distribution by:

Input = Numeric.arange(-1024,1024)/20.0
Output = Numeric.exp(-(Input**2))

At the edges of the distribution the largish negative arguments cause
the f.p. exception and prevent calculation of the entire array.

On windows the underflow errors are presumably rounded down to 0.0
automatically. Is it possible to get this behaviour on linux/ix86?

Anyone got any suggestions on a workarround for this? I see the numpy
docs suggest using Masked Arrays, but since I don't know in advance
which array items will generate the exceptions, I don't know which items
to mask. Can the MA module generate an array mask based on f.p.
exceptions from a calculation?

any suggestions would be much appreciated.

Bryan


-- 
Bryan Cole
Teraview Ltd., 302-304 Cambridge Science Park, Milton Road, Cambridge
CB4 0WG, United Kingdom.
tel: +44 (1223) 435380 / 435386 (direct-dial) fax: +44 (1223) 435382

Re: [Numpy-discussion] Newbie MA question

[Reggie D. <re...@me...>]

Cath,

The simple answer to your question is to use MA.set_print_limit(0) to
change the behavior of str and repr back to the standard Numeric
behavior of always printing the entire array.  Then you will get the
whole array, with the masked values as '--'.

Reggie

> Hi,
> 
> I'm not sure this is the right list - looking at the archives you guys 
> look more like developers more than users.
> 
> Any pointers to the proper place, or TFMs to R, are welcome. (Yes, I 
> have & have read the Numerical Python documentation.)
> 
> That said, I'm trying to use masked arrays (for a bioinformatics app 
> calculating stats on distances between pairs of amino acids, if you want 
> to know) and am having a little trouble printing.
> 
> One possible bug - "x.__str__()" and "x.__repr__()" show different 
> typecodes.
> Example:
>  >>> import MA
>  >>> id = MA.identity(20)
>  >>> r  = MA.ones((20,20))*1.0
>  >>> rm = MA.MaskedArray(r, mask=id)
>  >>> rm		       # ie, call x.__repr__()
> array(data =
>   array (20,20) , type = d, has 400 elements,
>        mask =
>   array (20,20) , type = 1, has 400 elements,
>        fill_value=[  1.00000002e+20,])
> 
>  >>> print rm         # ie, call x.__str__()
> array (20,20) , type = O, has 400 elements
>  >>>
> 
> Secondly, I can't see how to make this print the whole thing if I want 
> it. Which I do; all my matrices are 20*20 'cause there are 20 amino 
> acids.
> 
> The Numeric.array2string method can't be used because the data is always 
> masked.
> The best I've found so far is to make a Numeric array by filling with an 
> impossible number like -1, then use array2string. But visually it would 
> be nicer to see the masked values as "--".
> 
> Thanks for any help,
> regards,
> Cath
> Cath Lawrence,                       Cat...@an...
> Scientific Programmer,  Centre for Bioinformation Science,
> John Curtin School of Medical Research
> Australian National University,  Canberra ACT 0200
> ph: Maths (02) 6125 2904;  JCSMR (02) 6125 0417;
> mobile: 0421-902694   fax: (02) 61254712
> 
> 
> 
> -------------------------------------------------------
> This sf.net email is sponsored by:ThinkGeek
> Welcome to geek heaven.
> http://thinkgeek.com/sf
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion

[Numpy-discussion] Newbie MA question

[Cath L. <Cat...@an...>]

Hi,

I'm not sure this is the right list - looking at the archives you guys 
look more like developers more than users.

Any pointers to the proper place, or TFMs to R, are welcome. (Yes, I 
have & have read the Numerical Python documentation.)

That said, I'm trying to use masked arrays (for a bioinformatics app 
calculating stats on distances between pairs of amino acids, if you want 
to know) and am having a little trouble printing.

One possible bug - "x.__str__()" and "x.__repr__()" show different 
typecodes.
Example:
 >>> import MA
 >>> id = MA.identity(20)
 >>> r  = MA.ones((20,20))*1.0
 >>> rm = MA.MaskedArray(r, mask=id)
 >>> rm		       # ie, call x.__repr__()
array(data =
  array (20,20) , type = d, has 400 elements,
       mask =
  array (20,20) , type = 1, has 400 elements,
       fill_value=[  1.00000002e+20,])

 >>> print rm         # ie, call x.__str__()
array (20,20) , type = O, has 400 elements
 >>>

Secondly, I can't see how to make this print the whole thing if I want 
it. Which I do; all my matrices are 20*20 'cause there are 20 amino 
acids.

The Numeric.array2string method can't be used because the data is always 
masked.
The best I've found so far is to make a Numeric array by filling with an 
impossible number like -1, then use array2string. But visually it would 
be nicer to see the masked values as "--".

Thanks for any help,
regards,
Cath
Cath Lawrence,                       Cat...@an...
Scientific Programmer,  Centre for Bioinformation Science,
John Curtin School of Medical Research
Australian National University,  Canberra ACT 0200
ph: Maths (02) 6125 2904;  JCSMR (02) 6125 0417;
mobile: 0421-902694   fax: (02) 61254712

[Numpy-discussion] more functions for MLab.py: find, unique, union, intersect, etc.

[Mark H. <hi...@MI...>]

I've written some glue functions along the lines of what's in MLab.py:
 find
 unique
 union
 intersect
 setdiff
 setxor

Quick questions:
  * any opinions on the algorithms used?
  * are these worth adding to MLab.py?  If so, do I need to provide 
    documentation updates?  I can send this as a patch to MLab.py too.
  * Does anyone have opinions on what to do about the optional outputs from 
    matlab functions?  Is it worth trying to provide index matrices as second 
    arguments, for example?  If so, any opinions on implementation?

a matlab user trying to get started with Numeric,
        Mark

Re: [Numpy-discussion] Sparse matrices

[Jochen <jo...@un...>]

On Tue, 30 Jul 2002 09:42:13 -0400 Victor S Miller wrote:

Victor> I had noticed that Travis Oliphant had a sparse.py package,
Victor> but it no longer is available (clicking on the link gives a
Victor> "404").

It's part of scipy now.

Greetings,
Jochen
-- 
University of North Carolina                       phone: +1-919-962-4403
Department of Chemistry                            phone: +1-919-962-1579
Venable Hall CB#3290 (Kenan C148)                    fax: +1-919-843-6041
Chapel Hill, NC 27599, USA                            GnuPG key: 44BCCD8E

Re: [Numpy-discussion] Sparse matrices

[<vi...@id...>]

Sorry, I had a typo in the program.  It should be:

# M is n by k, and represents a sparse n by n matrix A
# the non-zero entries of row i of A start in column loc[i]
# and are the i-th row of M in locations loc[i]:loc[i]+k
# loc is the location vector
n,k = M.shape
mm = reshape(v,(-1,1))*M
w = zeros((n+m),v.typecode())
# is there a trick to replace the loop below?
for i in range(mm.shape[0]):
    w[loc[i]:loc[i]+k] += mm[i]
w = w[:n]

-- 
Victor S. Miller     | " ... Meanwhile, those of us who can compute can hardly
vi...@id...    | be expected to keep writing papers saying 'I can do the
CCR, Princeton, NJ   | following useless calculation in 2 seconds', and indeed
    08540 USA        | what editor would publish them?"  -- Oliver Atkin

[Numpy-discussion] Sparse matrices

[<vi...@id...>]

I had noticed that Travis Oliphant had a sparse.py package, but it no
longer is available (clicking on the link gives a "404").  I have a
particular kind of sparse matrix that I'd like to use to give vector
matrix multiplies.  In particular, it's an n x n matrix which has at
most k (which is small, usually 2 or 3) non-zeros in each row which
are in consecutive locations.  I have this encoded as an n x k matrix,
the i-th row gives the non-zero values in the i-th row of the big
matrix, and an n long vector of indices -- the i-th element gives the
starting position in the i-th row.  When I want to multiply this
matrix by a row vector v on the left.  To do the multiplication I do
the following:

# loc is the location vector
n = matrix.shape[0]
mm = reshape(v,(-1,1))*matrix
w = zeros((n+m),v.typecode())
for i in range(mm.shape[0]):
    w[loc[i]:loc[i]+matrix.shape[1]] += w[i]
w = w[:n]


I would like to be able to replace the loop with some Numeric
operations.  Is there a trick to do this?  Note that the n that I'm
using is around 100000, so that storing the full matrix is out of the
question (and multiplying by that matrix would be extremely
inefficient, anyway).
-- 
Victor S. Miller     | " ... Meanwhile, those of us who can compute can hardly
vi...@id...    | be expected to keep writing papers saying 'I can do the
CCR, Princeton, NJ   | following useless calculation in 2 seconds', and indeed
    08540 USA        | what editor would publish them?"  -- Oliver Atkin

[Numpy-discussion] A probable solution...

[Magnus L. H. <ma...@he...>]

After posting to the list (sorry about that ;) a possible solution
occurred to me... To get an approximation, I used fft(seq, 10) and
then inverted that using inverse_fft(signature, 100)... I guess that
fouled up the scale of things -- when I use fft(seq, 100)[:10] to get
the signature, it seems that everything works just fine...

Even though this _seems_ to do the right thing, I just wanted to make
sure that I'm not doing something weird here...

-- 
Magnus Lie Hetland                                  The Anygui Project
http://hetland.org                                  http://anygui.org

[Numpy-discussion] Spectral approximation/DFT

[Magnus L. H. <ma...@he...>]

Hi!

Sorry to ask what is probably a really clueless question -- if there
are any obvious sources of information about this, I'd be happy to go
there and find this out for myself... :]

Anyway; I'm trying to produce a graph to illustrate a time sequence
indexing method, which relies on extracting the low-frequent Fourier
coefficients and indexing a vector consisting of those. The graph
should contain the original time sequence, and one reconstructed from
the Fourier coefficients. Since it is reconstructed from only the
low-frequent coefficients (perhaps 10-20 coefficients), it will look
wavy and sinus'y.

Now... I'm no expert in signal processing (or the specifics of FFT/DFT
etc.), and I can't seem to make the FFT module do exactly what I want
here...

It seems that using fft(seq).real extracts the coefficients I'm after
(though I'm not sure whether the imaginary components ought to figure
in the equation somehow...)

But no matter how I use inverse_fft or inverse_real_fft it seems I
have to supply a number of coefficients equal to the sequence I want
to approximate -- otherwise there will be a huge offset between them.
Why is this so? Shouldn't the first coefficient take care of such an
offset? Perhaps inverse_fft isn't doing what I think it is?

If I haven't expressed myself clearly, I'd be happy to elaborate...

(For those who might be interested, the approach is described in the
paper found at http://citeseer.nj.nec.com/307308.html with a figure of
the type I'm trying to produce at page 5.)

Anyway, thanks for any help :)

-- 
Magnus Lie Hetland                                  The Anygui Project
http://hetland.org                                  http://anygui.org

Re: [Numpy-discussion] index method for array objects?

[Todd M. <jm...@st...>]

Aureli Soria Frisch wrote:

> Hi all,
>
> Has someone implemented a function for arrays that behaves like the 
> index(*) method for lists (it should then consider something like a 
> tolerance parameter).
>
> I suppose it could be maybe done with array.tolist() and list.index(), 
> but have someone implemented something more elegant/array-based?
>
> Thanks in advance
>
> Aureli
>
> PD: (*) index receive a value as an argument and retunrs the index of 
> the list member equal to this value...

I think the basics of what you're looking for are something like:

def index(a, b, eps):
      return nonzero(abs(a-b) < eps)

which should return all indices at which the absolute value of the 
difference between elements of a and b differ by less than eps.   e.g.:

 >>> import Numeric
 >>> index(Numeric.arange(10,20), 15, 1e-5)
array([5])

Todd

-- 
Todd Miller 			jm...@st...
STSCI / SSG			(410) 338 4576

[Numpy-discussion] index method for array objects?

[Aureli S. F. <au...@ip...>]

Hi all,

Has someone implemented a function for arrays that behaves like the 
index(*) method for lists (it should then consider something like a 
tolerance parameter).

I suppose it could be maybe done with array.tolist() and 
list.index(), but have someone implemented something more 
elegant/array-based?

Thanks in advance

Aureli

PD: (*) index receive a value as an argument and retunrs the index of 
the list member equal to this value...
-- 
#################################
Aureli Soria Frisch
Fraunhofer IPK
Dept. Pattern Recognition

post: Pascalstr. 8-9, 10587 Berlin, Germany
e-mail: au...@ip...
fon: +49 30 39006-143
fax: +49 30 3917517
web: http://vision.fhg.de/~aureli/web-aureli_en.html
#################################

Re: [Numpy-discussion] Numarray: Summary (seeding): personal code and manual suggestions on initial seeding in module RNG and RandomArray(2)

[Eric M. <e.m...@po...>]

Dear crunchers,

Please see my personal thoughts on the past discussion about initial
seeds some paragraphs down below, where I'd like to list concrete
code and manual enhancements aimed at providing users with a clear
understanding of it's usage (and pitfalls e.g. w/r to cryptographic
applications)...

=3D=3D> Suggestions for code and manual changes w/r to initial seeding
    (down below)

But first a response to Paul's earlier message:

On Thursday 25 July 2002 08:08, Paul F Dubois wrote:
> I'm not going to change the default seed on RNG. Existing users have th=
e
> right to stability, and not to have things change because someone think=
s
> a certain choice among several reasonable ones is better than the one
> previously made.

Well, I wasn't aware of the fact that things were completely set in stone
for Numarray solely for backward compatibilty. It was my impression that
numarray and it's accompanying xx2 packages were also open for redesign.
I agree stability is important, but numarray already breaks with Numeric
in other aspects so why should RNG (RNG2 in numarray?) or other packages
not be? It's more a matter of well documenting changes I think.
Users switching to numarray will already have to take into account some
changes and verify their code.
It's not that I "think a certain choice among several reasonable ones is
better" [although my favorite is still a fixed seed, as in RNG, for reaso=
ns
of reproducibility in later re-runs of Monte Carlo's that are not possibl=
e
now, because the naive user, using a clock seed, may not have saved the
initial seed with get_seed], but that the different packages, i.c. RNG
(RNG2 to be?) and RandomArray2, should be orthogonal in this respect.
I.e. the same, so 'default always an automagical (clock whatever) random
initial seed _or_ a fixed one'. Orthogonality is a very common and accept=
ed
design principle in computing science and for good reasons (usability).
Users changing from one PRNG to another (and using the default seed) woul=
d
otherwise be unwelcomely surprised by a sudden change in behavior of thei=
r
program.
I try to give logical arguments and real code examples in this discussion
and fail to see in Paul's reaction where I'm wrong.
By the way: in Python 2.1 alpha 2 seeding changed, too:
"""
- random.py's seed() function is new.  For bit-for-bit compatibility with
  prior releases, use the whseed function instead.  The new seed function
  addresses two problems:  (1) The old function couldn't produce more tha=
n
  about 2**24 distinct internal states; the new one about 2**45 (the best
  that can be done in the Wichmann-Hill generator).  (2) The old function
  sometimes produced identical internal states when passed distinct
  integers, and there was no simple way to predict when that would happen=
;
  the new one guarantees to produce distinct internal states for all
  arguments in [0, 27814431486576L).
"""

> There is the further issue here of RNG being advertised as similar to
> Cray's ranf() and that similarity extends to this default. Not to
> mention that for many purposes the current default is quite useful.

Perhaps I'm mistaken here, but RNG/Lib/__init__.py does
(-1 -> uses fixed internal seed):

    standard_generator =3D CreateGenerator(-1)

and:
    def ranf():
            "ranf() =3D a random number from the standard generator."
            return standard_generator.ranf()

And indeed Mixranf in RNG/Src/ranf.c does set them to 0:

    ...
    if(*s < 0){ /* Set default initial value */
        s48[0] =3D s48[1] =3D 0;
        Setranf(s48);
        Getranf(s48);

And this code, or I'm missing the point, uses a standard generator
from RNG, which demonstrates the same sequence of initial seeds in
re-runs (note that it does not suffer from the "1-second problem"
as RandomArray2 does, see the Appendix below for a demonstration
of that, because RNG uses milliseconds).
Note that 'ranf()' is listed in chapter 18 in Module RNG as one of
the 'Generator objects':

    $ python
    Python 2.2.1 (#1, Jun 25 2002, 20:45:02) ...
    >>> from numarray import *
    >>> from RNG import *
    >>> for i in range(3):
    ...    standard_generator.ranf()
    ...
    0.58011364857958725
    0.95051273498076583
    0.78637142533060356
    >>>

    $ python
    Python 2.2.1 (#1, Jun 25 2002, 20:45:02) ...
    >>> from numarray import *
    >>> from RNG import *
    >>> for i in range(3):
    ...     standard_generator.ranf()
    ...
    0.58011364857958725
    0.95051273498076583
    0.78637142533060356
    >>>


Ok, now then my own (and possibly biased) personal summary of the
past discussions and concrete code and manual recommendations:

=3D=3D> Suggestions for code and manual changes w/r to initial seeding

Conclusions:
1. Default initial seeding should be random (and not fixed).
   This is the general consensus and while it may not win the beauty
   contest in purist software engineering circles, it also is the
   default behavior in Python's own Random/WHRandom modules. URL:
       http://web.pydoc.org/2.2/random.html
=3D> Recommendations:
   - Like Python's random/whrandom module, default arguments to seed()
     should not be 0, but None, and this triggers the default behavior
     which is to use a random initial seed (ideally 'truly' random from
     e.g. /dev/random or otherwise clock or whatever based), because:
     o better usability: users changing from Python's own random to
       numarray's random facilities will find familiar seed() usage
       semantics
     o often 0 itself can be a legal seed (although the MersenneTwister
       does not recommend it)
   - Like RNG provide support for using a built-in fixed seed by
     supplying negative seeds to seed(), rationale:
     o support for reproducible re-runs of Monte Carlo's without
       having to specify ones own initial seed
     o usability: naive users may not know a 'good' seed is, like:
       can it be 0 or must it be >0, what is the maximum, etc.
   - See my suggested code fix for RandomArray2.seed() in the Appendix be=
low.
   - Likewise, in RNG:
     o CreateGenerator (s, ...) should be changed to CreateGenerator (s=3D=
None)
       Also note Python's own: def create_generators(num, delta, firstsee=
d=3DNone)
       from random (random.py), url: http://web.pydoc.org/2.2/random.html
     o RNG's code should be changed from testing on 0 to testing on None =
first
       (which results in using the clock), then on < 0 (use built-in seed=
),
       and then using the user provided seed (which is thus >=3D 0, and h=
ence
       can also be 0)
     o 'standard_generator =3D CreateGenerator(-1)' should be changed to
       'standard_generator =3D CreateGenerator() and results in using the=
 clock
   - Put some explicit warnings in the numarray manual, that the seeding
     of numarray's packages should _not_ be used in those parts of softwa=
re
     where unpredictability of seeds is important, such as for example,
     cryptographical software for creating session keys, TCP sequence num=
bers
     etc. Attacks on crypto software usually center around these issues.
     Ideally, a /dev/random should be used, but with the current system
     clock based implementation, the seeds are not random, because the cl=
ock
     does not have deci-nanosecond precision (10**10 ~=3D 2**32) yet ;-)

Appendix
--------
** 1. "1-second problem" with RandomArray2:

$ python
Python 2.2.1 (#1, Jun 25 2002, 20:45:02) ...
>>> from numarray import *
>>> from RandomArray2 import *
>>> import time
>>> import sys
>>> sys.version
'2.2.1 (#1, Jun 25 2002, 20:45:02) \n[GCC 2.95.3 20010315 (SuSE)]'
>>> numarray.__version__
'0.3.5'
>>> for i in range(3):
=2E..     time.time()
=2E..     RandomArray2.seed()
=2E..     RandomArray2.get_seed()
=2E..     time.sleep(1)
=2E..     print
=2E..
1027591910.9043469
(102759, 1911)

1027591911.901091
(102759, 1912)

1027591912.901088
(102759, 1913)

>>> for i in range(3):
=2E..     time.time()
=2E..     RandomArray2.seed()
=2E..     RandomArray2.get_seed()
=2E..     time.sleep(0.3)
=2E..     print
=2E..
1027591966.260392
(102759, 1967)

1027591966.5510809
(102759, 1967)

1027591966.851079
(102759, 1967)

Note that Python (at least 2.2.1) own random() suffers much less
from this (on my 450 MHz machine, every 10-th millisecond or so
the seed will be different):

$ python
Python 2.2.1 (#1, Jun 25 2002, 20:45:02) ...
>>> from random import *
>>> import time
>>>
>>> for i in range(3):
=2E..     print long(time.time() * 256)
=2E..
263065231349
263065231349
263065231349
>>> for i in range(3):
=2E..     print long(time.time() * 256)
=2E..     time.sleep(.00001)
=2E..
263065240314
263065240315
263065240317

By the way, Python's own random.seed() also suffers from this,
but on a 10th-millisecond level (on my 450 Mhz i586 at least).
For the implementation of seed() see Lib/random.py, basically a
'long(time.time()' is used:

$ python
Python 2.2.1 (#1, Jun 25 2002, 20:45:02) ...
>>> from random import *
>>> import time
>>> for i in range(3):
=2E..     print long(time.time() * 256)
=2E..
263065231349
263065231349
263065231349
>>> for i in range(3):
=2E..     print long(time.time() * 256)
=2E..     time.sleep(.00001)
=2E..
263065240314
263065240315
263065240317

2. Proposed re-implementation of RandomArray2.seed():

def seed(x=3DNone,y=3DNone):
    """seed(x, y), set the seed using the integers x, y:
       x or y is None (or not specified):
          A random seed is used which in the current implementation
          may be based on the system's clock.
          Warning: do not this seed in software where the initial seed ma=
y
          not be predictable, such as for example, in cryptographical sof=
tware
          for creating session keys.
       x < 0 or y < 0:
          Use the module's fixed built-in seed which is the tuple
          (1234567890L, 123456789L) (or whatever)
       x >=3D 0 and y >=3D 0
          Use the seeds specified by the user.
          (Note: some random number generators do not recommend using 0)
       Note: based on Python 2.2.1's random.seed(a=3DNone).
       ADAPTED for _2_ seeds as required by ranlib.set_seeds(x,y)
    """
    if (x !=3D None and type (x) !=3D IntType) or
       (y !=3D None and type (y) !=3D IntType) :
        raise ArgumentError, "seed requires integer arguments (or None)."
    if x =3D=3D None or y =3D=3D None:
        try:
            # This would be the best, but is problematic under Windows/Ma=
c.
            # To my knowledge there isn't a portable lib_randdevice yet.
            # As GPG, OpenSSH and OpenSSL's code show, getting entropy
            # under Windows is problematic.
            # However, Python 2.2.1's socketmodule does wrap the ssl code=
=2E
            import dev_random_device  # uses /dev/random or equivalent
            x =3D dev_random_device.nextvalue()   # egd.sf.net is a user =
space
            y =3D dev_random_device.nextvalue()   # alternative
        except:
            # Use Mixranf() from RNG/Src/ranf.c here or, perhaps better,
            # use Python 2.2.1's code? At least it looks simpler and does=
 not
            # have the platform dependency's and has possibly met wider t=
esting
            # (and why not re-use code? ;-)
            # For Python 2.2.1's random.seed(a=3DNone), see url:
            #     http://web.pydoc.org/2.2/random.html
            # and file Lib/random.py.
            # Do note, however, that on my 450 Mhz machine, the statement
            # 'long(time.time() * 256)' will generate the same values
            # within a tenth of a millisecond (see Appendix #1 for a code
            # example). This can be fixed by doing a time.sleep(0.001).
            # See my #EM# comment.
            # Naturally this code needs to be adapted for ranlib's
            # generator, because this code uses the Wichmann-Hill generat=
or.

---cut: Wichmann-Hill---
    def seed(self, a=3DNone):
        """Initialize internal state from hashable object.

        None or no argument seeds from current time.

        If a is not None or an int or long, hash(a) is used instead.

        If a is an int or long, a is used directly.  Distinct values betw=
een
        0 and 27814431486575L inclusive are guaranteed to yield distinct
        internal states (this guarantee is specific to the default
        Wichmann-Hill generator).
        """

        if a is None:
            # Initialize from current time
            import time
            a =3D long(time.time() * 256)
            #EM# Guarantee unique a's between subsequent call's of seed()
            #EM# by sleeping one millisecond. This should not be harmful,
            #EM# because ordinarily, seed() will only be called once or s=
o
            #EM# in a program.
            time.sleep(0.001)

        if type(a) not in (type(3), type(3L)):
            a =3D hash(a)

        a, x =3D divmod(a, 30268)
        a, y =3D divmod(a, 30306)
        a, z =3D divmod(a, 30322)
        self._seed =3D int(x)+1, int(y)+1, int(z)+1
---cut: Wichmann-Hill---
    elif x < 0 or y < 0:
        x =3D 1234567890L  # or any other suitable 0 - 2**32-1
        y =3D 123456789L
    ranlib.set_seeds(x,y)

3. Mersenne Twister, another PRNG:



Bye-bye,

Eric
--=20
Eric Maryniak <e.m...@po...>
WWW homepage: http://pobox.com/~e.maryniak/
Mobile phone: +31 6 52047532, or (06) 520 475 32 in NL.

In a grocery store, the Real Programmer is the one who insists on running
the cans past the laser checkout scanner himself, because he never could
trust keypunch operators to get it right the first time.

RE: [Numpy-discussion] Numarray: question on RandomArray2.seed(x=0, y=0) system clock default and possible bug

[Paul F D. <pa...@pf...>]

I'm not going to change the default seed on RNG. Existing users have the
right to stability, and not to have things change because someone thinks
a certain choice among several reasonable ones is better than the one
previously made. 

There is the further issue here of RNG being advertised as similar to
Cray's ranf() and that similarity extends to this default. Not to
mention that for many purposes the current default is quite useful.

Re: [Numpy-discussion] Numarray: question on RandomArray2.seed(x=0, y=0) system clock default and possible bug

[<pet...@no...>]

Just to stick my oar in:

I think Eric's preference is predicated by the lousiness (or otherwise?)  
of RandomArray's seeding mechanism. The random sequences generated by
incremental seeds should, by design, be uncorrelated thus allowing the use
of the system clock as a seed source.

If you're running lots of simulations (as I do with Monte Carlos, though
not in numpy) using PRNGs, the last thing you want is the task to find a
(pseudo) random source of seeds. Using /dev/random is not particularly
portable; the system clock is much easier to obtain and is fine as long as
your iteration cycle is longer than its resolution.

Peter