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

Re: [Numpy-discussion] Bug: extremely misleading array behavior

[eric <er...@en...>]

Hey Larry,

I actually thought, as you did, that indexing the array returns an element
converted to a scalar -- and it does in the "default" cases when you don't
specify a non-standard typecode.

After testing, it looks like values that are representable as native Python
types ('l', 'd', and 'D') are returned as actual values while non-standard types
are returned as views into the array.  Is this intentional?  It is dangerous to
have the behavior change based on the type.  It seems they should all be views
or they should all  be converted to a scalar.

Here is your test code modified to test all Numeric types:

    import Numeric

    def test_index(typecode):
        print 'typcode:', typecode
        a = Numeric.zeros((2, 2), typecode)
        n = a[1, 1]     # fetch interesting value from array
        print n
        a[1, 1] = 10    # change array
        print n         # blam
        print type(n)   # huh
        print

    print 'Numeric version:', Numeric.__version__
    for t in ['i','1','s','l','f','d','F','D']:
        test_index(t)

And here is the output.  Look at the types returned.

C:\home\ej\wrk\junk>python num_index.py
Numeric version: 21.0
typcode: i
0
10
<type 'array'>

typcode: 1
0
10
<type 'array'>

typcode: s
0
10
<type 'array'>

typcode: l
0
0
<type 'int'>

typcode: f
0.0
10.0
<type 'array'>

typcode: d
0.0
0.0
<type 'float'>

typcode: F
0j
(10+0j)
<type 'array'>

typcode: D
0j
0j
<type 'complex'>

eric

----- Original Message -----
From: "Larry Denneau" <la...@pa...>
To: <num...@li...>
Sent: Tuesday, May 28, 2002 1:13 PM
Subject: [Numpy-discussion] Bug: extremely misleading array behavior


> Hello,
>
> I recently discovered the following behavior when fetching values
> from a Numeric array.  Can somebody offer some insight?
>
> #1)
>
>   import Numeric
>
>   a = Numeric.zeros((2, 2), 'i')
>   n = a[1, 1]     # fetch interesting value from array
>   print n
>   a[1, 1] = 10    # change array
>   print n         # blam
>   print type(n)   # huh
>
>   [bash]$ python 1.py
>   0
>   10
>   <type 'array'>
>
> but
>
> #2)
>
>   import Numeric
>
>   a = Numeric.zeros((2,), 'i')
>   n = a[1]
>   print n
>   a[1] = 10
>   print n
>   print type(n)
>
>   [bash]$ python 2.py
>   0
>   0
>   <type 'int'>
>
> #2 works the way one would expect, and #1 does not (n changes).
> They should at least both behave the same. :-)  At a minimum, naive
> use of arrays can lead to confusing or disastrous results, since
> a single value fetched from an array can change behind your back.
>
> It appears n is aliased into a, but preserves its value when a is
> deleted (with del(a)).  What happens to the "rest of" a?
>
> I'm using Python 2.2, Numeric-21.0, on both Unix and Win32.
>
> Thanks,
> Larry
>
> _______________________________________________________________
>
> Don't miss the 2002 Sprint PCS Application Developer's Conference
> August 25-28 in Las Vegas -- http://devcon.sprintpcs.com/adp/index.cfm
>
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>

Re: [Numpy-discussion] Bug: extremely misleading array behavior

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

Hi Larry,

On Tue, 28 May 2002, Larry Denneau wrote:

> All the Numpy documentation examples (see
> http://pfdubois.com/numpy/html2/numpy-6.html#pgfId-36033, "Getting and
> Stting Array Values") use the [x, y] notation instead of [x][y], so I
> would consider this a bug in the documentation, since the [x, y]
> method leads to unexpected behavior.

If you look the section "Slicing Arrays" then a[1] is actually a[1,:],
that is, an one dimensional array. From your description, a[1,1] must be
an array with 0 rank. It seems that the Numeric documentation is missing
(though, I didn't look too hard) the following rules of thumb:

  If `a' is rank 1 array, then a[i] is Python scalar or object. [MISSING]
  If `a' is rank > 1 array, then a[i] is a sub-array a[i,...]

> I'm still curious what happens to the original array when 
> 
>   n=a[1, 1]
>   del(a)

I think the original array `a' is not actually deleted until `n' gets
deleted. If I recall correctly, then `n' is a sub-array of `a' so that
internally it contains only a reference to `a' in the sense that
a.data==n.data but strides and dimension arrays differ.

Pearu

Re: [Numpy-discussion] Bug: extremely misleading array behavior

[Larry D. <la...@pa...>]

Pearu Peterson said:
> 
> On Tue, 28 May 2002, Larry Denneau wrote:
> 
> > Hello,
> > 
> > I recently discovered the following behavior when fetching values
> > from a Numeric array.  Can somebody offer some insight?
> > 
> > #1)
> > 
> >   import Numeric
> > 
> >   a = Numeric.zeros((2, 2), 'i')
> >   n = a[1, 1]     # fetch interesting value from array
> >   print n 
> >   a[1, 1] = 10    # change array
> >   print n         # blam
> >   print type(n)   # huh
> > 
> >   [bash]$ python 1.py
> >   0
> >   10
> >   <type 'array'>

[ deleted]

> Use
>      
>     a[1][1] = 10
>      
> and the output will be
>   
> 0
> 0 
> <type 'int'>
> 
> I find it is an useful feature in Numeric to have both behaviours of
> either using a[1,1] or a[1][1]. You may want to dig into Numeric's
> userguide to get a more detailed explanation of the differences.
> 
> Regards,
> 	Pearu

Hi Pearu,

I assume you mean

  n = a[1][1]

which produces the expected behavior.  All the Numpy documentation
examples (see http://pfdubois.com/numpy/html2/numpy-6.html#pgfId-36033,
"Getting and Stting Array Values") use the [x, y] notation instead
of [x][y], so I would consider this a bug in the documentation,
since the [x, y] method leads to unexpected behavior.

I'm still curious what happens to the original array when 

  n=a[1, 1]
  del(a)

but that may have to wait until I have time to peruse the Numeric
source.

Thanks,
Larry

Re: [Numpy-discussion] Bug: extremely misleading array behavior

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

On Tue, 28 May 2002, Pearu Peterson wrote:

<snip>

> Use
>      
>     a[1][1] = 10
>      

Oops, I meant

  n = a[1][1]

Pearu

Re: [Numpy-discussion] Bug: extremely misleading array behavior

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

On Tue, 28 May 2002, Larry Denneau wrote:

> Hello,
> 
> I recently discovered the following behavior when fetching values
> from a Numeric array.  Can somebody offer some insight?
> 
> #1)
> 
>   import Numeric
> 
>   a = Numeric.zeros((2, 2), 'i')
>   n = a[1, 1]     # fetch interesting value from array
>   print n 
>   a[1, 1] = 10    # change array
>   print n         # blam
>   print type(n)   # huh
> 
>   [bash]$ python 1.py
>   0
>   10
>   <type 'array'>
> 
> but
> 
> #2)
> 
>   import Numeric
> 
>   a = Numeric.zeros((2,), 'i')
>   n = a[1]
>   print n
>   a[1] = 10
>   print n
>   print type(n)
> 
>   [bash]$ python 2.py
>   0
>   0
>   <type 'int'>
> 
> #2 works the way one would expect, and #1 does not (n changes).
> They should at least both behave the same. :-)  At a minimum, naive
> use of arrays can lead to confusing or disastrous results, since
> a single value fetched from an array can change behind your back.

Use
     
    a[1][1] = 10
     
and the output will be
  
0
0 
<type 'int'>

I find it is an useful feature in Numeric to have both behaviours of
either using a[1,1] or a[1][1]. You may want to dig into Numeric's
userguide to get a more detailed explanation of the differences.

Regards,
	Pearu

[Numpy-discussion] Bug: extremely misleading array behavior

[Larry D. <la...@pa...>]

Hello,

I recently discovered the following behavior when fetching values
from a Numeric array.  Can somebody offer some insight?

#1)

  import Numeric

  a = Numeric.zeros((2, 2), 'i')
  n = a[1, 1]     # fetch interesting value from array
  print n 
  a[1, 1] = 10    # change array
  print n         # blam
  print type(n)   # huh

  [bash]$ python 1.py
  0
  10
  <type 'array'>

but

#2)

  import Numeric

  a = Numeric.zeros((2,), 'i')
  n = a[1]
  print n
  a[1] = 10
  print n
  print type(n)

  [bash]$ python 2.py
  0
  0
  <type 'int'>

#2 works the way one would expect, and #1 does not (n changes).
They should at least both behave the same. :-)  At a minimum, naive
use of arrays can lead to confusing or disastrous results, since
a single value fetched from an array can change behind your back.

It appears n is aliased into a, but preserves its value when a is
deleted (with del(a)).  What happens to the "rest of" a?

I'm using Python 2.2, Numeric-21.0, on both Unix and Win32.

Thanks,
Larry

Re: [Numpy-discussion] assignment to matrix slice in Numeric bug/feature?

[Roman S. <rn...@on...>]

On Tue, 21 May 2002, Roman Suzi wrote:

>hello,
>
>I've found that the following fragment of code gives an error while with
>other shapes of b there is no problem:

...

> The mentioned behaviour looks like a bug, because
> a[1:3,1:2] and b have the same shape and according to docs
> b must be copied to a-slice 1:1...
>
> Numeric is version 21.0, Python 2.1 under Linux RedHat 7.2.
>
> Thank you in advance!

For the record:

- the bug is fixed in CVS.

Thank you, guys!

Sincerely yours, Roman Suzi
-- 
\_ Russia \_ Karelia \_ Petrozavodsk \_ rn...@on... \_
\_ Saturday, May 25, 2002 \_ Powered by Linux RedHat 7.2 \_
\_ "... All the world's a stage, and I missed rehearsal." \_

Re: [Numpy-discussion] Bug in MLab.std()

[John J. L. <jj...@po...>]

On Thu, 23 May 2002, Jonathan M. Gilligan wrote:

> I submitted this to the Sourceforge bug tracker, but wanted also to let this
> list know, as this is a potentially nasty bug.
>
> MLab.std() gives completely incorrect answers for multidimensional arrays
> when axis != 0.
[...]

There was a earlier bug with this function -- is there a regression test
to make sure this change doesn't fail the way it did before?


John

Re: [Numpy-discussion] Installation woes

[Paul B. <Ba...@st...>]

Eric Hagemann wrote:

> Hey All,
> 
>     Doing my first install on a linux system (Mandrake 8.2 on a PPC).
> 
> I downloaded the *.tar.gz and unzip/untarred it.  I enter the command
> 'python setup.py install' and get the following error
> 
> "open '/usr/lib/python2.2/config/Makefile' (No such file or directory)"
> 
> This is the python that comes with mandrake.  I have been able to download
> and install a new version of python2.2 and sucessfully install Numpy there
> but I was wondering what I needed to do to get it running in the 'base'
> python install
> 
> any pointers ?

Yes, you need to install the libpython2.2-devel RPM.  It contains this file.

-- 
Paul Barrett, PhD      Space Telescope Science Institute
Phone: 410-338-4475    ESS/Science Software Group
FAX:   410-338-4767    Baltimore, MD 21218

[Numpy-discussion] Installation woes

[Eric H. <eha...@co...>]

Hey All,

    Doing my first install on a linux system (Mandrake 8.2 on a PPC).

I downloaded the *.tar.gz and unzip/untarred it.  I enter the command
'python setup.py install' and get the following error

"open '/usr/lib/python2.2/config/Makefile' (No such file or directory)"

This is the python that comes with mandrake.  I have been able to download
and install a new version of python2.2 and sucessfully install Numpy there
but I was wondering what I needed to do to get it running in the 'base'
python install

any pointers ?

Cheers
Eric

[Numpy-discussion] Bug in MLab.std()

[Jonathan M. G. <jon...@va...>]

I submitted this to the Sourceforge bug tracker, but wanted also to let this
list know, as this is a potentially nasty bug.

MLab.std() gives completely incorrect answers for multidimensional arrays
when axis != 0.

>>> foo
array([[[ 1.,  1.,  1.],
        [ 2.,  2.,  2.],
        [ 3.,  3.,  3.]],
       [[ 1.,  4.,  4.],
        [ 2.,  5.,  5.],
        [ 3.,  6.,  6.]]])
>>> std(foo)
array([[ 0.        ,  2.12132034,  2.12132034],
       [ 0.        ,  2.12132034,  2.12132034],
       [ 0.        ,  2.12132034,  2.12132034]])
>>> std(foo, 1)
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])

The following should fix the problem (but I haven't tested it extensively):

def std(m,axis=0):
    """std(m,axis=0) returns the standard deviation along the given
    dimension of m.  The result is unbiased with division by N-1.
    If m is of integer type returns a floating point answer.
    """
    x = asarray(m)
    n = float(x.shape[axis])
    x2 = mean(x * x, axis)
    x = mean(x, axis)
    return sqrt((x2 - x * x) * n /(n-1.0))

Jonathan Gilligan

[Numpy-discussion] RE: Numpy-discussion digest, Vol 1 #463 - 1 msg

[Tariq R. <tar...@li...>]

 i'm interested in this too?

 in fact i am willing to help write some wavelet tools for python - but
you'll have to be patient as i've only just started on my PyObjects et al
... and i'm not an expert either - just very interested and entheusiastic!

 Tariq Rashid

-----Original Message-----
From: num...@li...
[mailto:num...@li...]On Behalf Of
num...@li...
Sent: 22 May 2002 20:05
To: num...@li...
Subject: Numpy-discussion digest, Vol 1 #463 - 1 msg


Send Numpy-discussion mailing list submissions to
	num...@li...

To subscribe or unsubscribe via the World Wide Web, visit
	https://lists.sourceforge.net/lists/listinfo/numpy-discussion
or, via email, send a message with subject or body 'help' to
	num...@li...

You can reach the person managing the list at
	num...@li...

When replying, please edit your Subject line so it is more specific
than "Re: Contents of Numpy-discussion digest..."


Today's Topics:

   1. Q: Wavelet tools with Python (Zaur Shiboukhov)

--__--__--

Message: 1
Date: Wed, 22 May 2002 15:27:49 +0400
From: Zaur Shiboukhov <sz...@fr...>
To: num...@li...
Organization: RI Applied Mathematics and Automation RAS
Subject: [Numpy-discussion] Q: Wavelet tools with Python

How is about wavelet tools with Python?

Zaur




--__--__--

_______________________________________________
Numpy-discussion mailing list
Num...@li...
https://lists.sourceforge.net/lists/listinfo/numpy-discussion


End of Numpy-discussion Digest

[Numpy-discussion] Q: Wavelet tools with Python

[Zaur S. <sz...@fr...>]

How is about wavelet tools with Python?

Zaur

[Numpy-discussion] assignment to matrix slice in Numeric bug/feature?

[Roman S. <rn...@on...>]

hello,

I've found that the following fragment of code gives an error while with
other shapes of b there is no problem:

----------------------------
#!/usr/bin/python2

from Numeric import *
from Matrix import *
a = Matrix(zeros([4,4]))
b = Matrix(ones([2,1]))
print a, a.shape
print b, b.shape
q = a[1:3,1:2]
print q, q.shape
a[1:3,1:2] = b
-----------------------------

resulting in:
-------------------------------------
Matrix([[0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0]]) (4, 4)
Matrix([[1],
       [1]]) (2, 1)
Matrix([[0],
       [0]]) (2, 1)
Traceback (most recent call last):
  File "./numpybug.py", line 11, in ?
    a[1:3,1:2] = b
  File "/usr/lib/python2.1/site-packages/Numeric/Matrix.py", line 180, in 
__setitem__
    def __setitem__(self, index, value): self.array[index] = 
asarray(squeeze(value),self._typecode)
ValueError: matrices are not aligned for copy
-----------------------------------------------

The mentioned behaviour looks like a bug, because
a[1:3,1:2] and b have the same shape and according to docs
b must be copied to a-slice 1:1...

Numeric is version 21.0, Python 2.1 under Linux RedHat 7.2.

Thank you in advance!

Sincerely yours, Roman Suzi
-- 
\_ Russia \_ Karelia \_ Petrozavodsk \_ rn...@on... \_
\_ Saturday, May 18, 2002 \_ Powered by Linux RedHat 7.2 \_

[Numpy-discussion] [ANN] Chaco: A 2D scientific graphics library

[eric <er...@en...>]

I'm happy to announce that Enthought is developing a platform independent
plotting library for Python. The Chaco project, as it is named, is funded by
the Space Telescope Science Institute (STScI) and licensed under a BSD style
open source license. Chaco is designed for presentation quality scientific 2D
graphics on a variety of output devices. The initial targets are wxPython,
TkInter, Mac OS X, and PDF for hard copy output. It's design is extensible so
that other backends, such as OpenGL, can be added. Currently, the low-level API
for wxPython, Mac OS X, and PDF are operational.  The high level graphics
objects will be developed over the coming months.  Chaco is hosted at the SciPy
site. For more information visit:

        http://www.scipy.org/site_content/chaco

People are invited to comment on and contribute to the project.  Chaco's
discussion
list is:

    sci...@sc...

To subscribe, go to the mailing list's info page:

    http://scipy.net/mailman/listinfo/scipy-chaco

thanks,
eric jones

--
Eric Jones <eric at enthought.com>
Enthought, Inc. [www.enthought.com and www.scipy.org]
(512) 536-1057

[Numpy-discussion] Some problems with Numeric documenation

[Edward C. J. <edc...@er...>]

Here are some documentation problems with Numeric (20.2.0). I use the html
form of the documentation.

1. descr->elsize is not documented.

2. The documentation for PyArray_FromDimsAndData says

       This function should only be used to access global data that will
       never be freed (like FORTRAN common blocks).

   The document writer copied this from the source code internal
   documentation.  As far as I can tell, this is wrong. I suggest

       The object returned by this function contains a pointer to a
       pre-existing block of data. If you delete the data before the
       Python reference count has dropped to zero, you will have a dangling
       pointer which is usually a disaster. One safe way to do this is to
       Py_DECREF the object in the same C function where it was created.

3. By a semi-automated application of grep and Python, I have found the
   following functions which appear to be in the C API and are not
   documented. Note that PyArray_FromDimsAndData and PyArray_CopyArray
   can be combined to do the C level equivalent of
       arr[a:b,c:d] = something

  PyArray_ArgMax
       multiarraymodule.c:extern PyObject *PyArray_ArgMax(PyObject *op) {
   PyArray_ArgSort
       multiarraymodule.c:extern PyObject *PyArray_ArgSort(PyObject *op) {
   PyArray_BinarySearch
       multiarraymodule.c:extern PyObject *PyArray_BinarySearch(PyObject *op1, PyObject *op2) {
   PyArray_CONTIGUOUS
       arrayobject.c:#define PyArray_CONTIGUOUS(m) (ISCONTIGUOUS(m) ? Py_INCREF(m), m : \
   PyArray_Choose
       multiarraymodule.c:extern PyObject *PyArray_Choose(PyObject *ip, PyObject *op) {
   PyArray_Concatenate
       multiarraymodule.c:extern PyObject *PyArray_Concatenate(PyObject *op) {
   PyArray_Converter
       arrayobject.c:extern int PyArray_Converter(PyObject *object, PyObject **address) {
   PyArray_CopyArray
       arrayobject.c:int PyArray_CopyArray(PyArrayObject *dest, PyArrayObject *src) {
   PyArray_CopyObject
       arrayobject.c:int PyArray_CopyObject(PyArrayObject *dest, PyObject *src_object) {
   PyArray_Correlate
       multiarraymodule.c:extern PyObject *PyArray_Correlate(PyObject *op1, PyObject *op2, int mode) {
   PyArray_FromDimsAndDataAndDescr
       arrayobject.c:PyObject *PyArray_FromDimsAndDataAndDescr(int nd, int *d, 
   PyArray_FromScalar
       arrayobject.c:PyObject *PyArray_FromScalar(PyObject *op, int type) {
   PyArray_InnerProduct
       multiarraymodule.c:extern PyObject *PyArray_InnerProduct(PyObject *op1, PyObject *op2) {
   PyArray_Item
       arrayobject.c:extern PyObject * PyArray_Item(PyObject *op, int i) {
   PyArray_NBYTES
       arrayobject.h:#define PyArray_NBYTES(mp) ((mp)->descr->elsize * PyArray_SIZE(mp))
   PyArray_Put
       arrayobject.c:extern PyObject *PyArray_Put(PyObject *self0, PyObject *indices0, 
   PyArray_PutMask
       arrayobject.c:extern PyObject *PyArray_PutMask(PyObject *self0, PyObject *mask0, 
   PyArray_Repeat
       multiarraymodule.c:extern PyObject *PyArray_Repeat(PyObject *aop, PyObject *op, int axis) {
   PyArray_Resize
       arrayobject.c:static PyObject * PyArray_Resize(PyArrayObject *self, PyObject *shape) {
   PyArray_Sort
       multiarraymodule.c:extern PyObject *PyArray_Sort(PyObject *op) {
   PyArray_TYPES
       arrayobject.h:enum PyArray_TYPES { PyArray_CHAR, PyArray_UBYTE, PyArray_SBYTE,
   PyArray_ToList
       arrayobject.c:static PyObject *PyArray_ToList(PyObject *self) {
   PyArray_Transpose
       multiarraymodule.c:extern PyObject *PyArray_Transpose(PyArrayObject *ap, PyObject *op) {
   PyArray_compare_lists
       arrayobject.c:extern int _PyArray_compare_lists(int *l1, int *l2, int n) {



       arrayobject.c:extern int _PyArray_compare_lists(int *l1, int *l2, int n) {

       multiarraymodule.c:extern PyObject *PyArray_Sort(PyObject *op) {

[Numpy-discussion] Pyfort version 8.0a1 in CVS

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

Pyfort version 8.0a1 is now in CVS for testing.

A new user interface helps scientists create packages. Packages can be
uninstalled. A GUI interface is provided for creating new packages and
managing their compilation and 
installation.

Documentation for the new features will be available shortly at
pyfortran.sf.net.

Re: [Numpy-discussion] Noise Simulation

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

On Tue, 14 May 2002, Johan Fredrik =D8hman wrote:

> Sure,
> Kasdin, N.J. and Walter, T. (1992), Discrete Simulation of Power Law No=
ise,
> 1992, IEEE Frequency Control Symposium, p. 274-283.
>=20
> This C program shown there, and thus my python program returns an finit=
e
> vector of noise.   Now, let's assume that I want more noise, how can I
> concatenate another vector ?   In the case of PWH (Phase white noise) t=
his
> is quite simple, but the other cases are more obscure.  Anybody knows w=
hat
> is a valid procedure to create "more noise"...

Look in "Fast Algorithms for DSP" by Richard E Blahut, p284
(Addison-Wesley, Mass. 1985) ISBN 0-201-10155-6 or "Numerical Recipes in
C", p543 (CUP, 2nd Ed, 1992).=20

What you're doing is convolving a random stream with a FIR filter via=20
FFTs. To apply to a larger stream, you need to convolve by sections using=
=20
either the overlap-save or the overlap-add method.

Regards,
 Peter

Re: [Numpy-discussion] Noise Simulation

[<jo...@oh...>]

Sure,
Kasdin, N.J. and Walter, T. (1992), Discrete Simulation of Power Law Noise,
1992, IEEE Frequency Control Symposium, p. 274-283.

This C program shown there, and thus my python program returns an finite
vector of noise.   Now, let's assume that I want more noise, how can I
concatenate another vector ?   In the case of PWH (Phase white noise) this
is quite simple, but the other cases are more obscure.  Anybody knows what
is a valid procedure to create "more noise"...

--
JFØ


----- Original Message -----
From: "Hartley, Ed" <e.h...@la...>
To: <num...@li...>
Sent: Tuesday, May 14, 2002 10:03 AM
Subject: [Numpy-discussion] Noise Simulation


> Hi
> interesting discussion please could you quote the complete reference to
> this
> > I have based my routine on the document "Discrete simulation of power
low
> > noise", IEEE International Frequency Control Symposium 27-29
> >
> Regards
> Ed Hartley
> Research Fellow
> Computing Department
> Lancaster University
> Lancaster
> UK LA1 4YR
> Phone +44 (0) 1524 593675
> Fax   +44 (0) 1524 593608
>
>
> _______________________________________________________________
>
> Have big pipes? SourceForge.net is looking for download mirrors. We supply
> the hardware. You get the recognition. Email Us: ban...@so...
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>

[Numpy-discussion] Noise Simulation

[Hartley, E. <e.h...@la...>]

	Hi 
	interesting discussion please could you quote the complete reference to
this 
> I have based my routine on the document "Discrete simulation of power low
> noise", IEEE International Frequency Control Symposium 27-29
> 
	Regards
	Ed Hartley
Research Fellow
Computing Department
Lancaster University
Lancaster 
UK LA1 4YR
Phone +44 (0) 1524 593675
Fax   +44 (0) 1524 593608

Re: [Numpy-discussion] real_fft

[<jo...@oh...>]

First, I would like to say I have solved the problem myself, but thanks to
those who tried to help me :).  The error was related to the way the NR
realft() function interpreted the parameters.



realft(data[],length,1)
I thought length was how long the data array is.  However the routine reads
and uses data
TWICE the length !! :-/



Now, back to the side note.  Yes, colored noise is quite interesting.  It
has applications especially in clock simulations.

I have based my routine on the document "Discrete simulation of power low
noise", IEEE International Frequency Control Symposium 27-29

This is a document not found on the internet, but if you ask your university
or library they might have it.   In the appendix of this document there is a
sample C program which simulates power law noise.


Now that my python port of this is working you can have the code :)
The parameters to the noiseGen are:
points = number of points to generate
X      = the array in which to add the noise¨
Qd    = The noise variance
b       = the power law variable.   f^(b + 2).  Where f is the frequency.
           i.e b = 0 gives what is called white phase noise
                b=-1  gives white flicker noise
                b=-2  white frequency noise (phase random walk noise)
                b=-3  flicker frequency noise (pink noise)
                b=-4  random walk frequency noise
Non-integer values of b is also allowed, .i.e -2.5

############################################################################
##
# MODULE NAME:   Colored Noise Module
# AUTHOR: Johan Fredrik Øhman
#
# This module should simulate the colored noise found in clocks
#
############################################################################
##
import math, FFT, RNG, sys, emath, Numeric

__gausian   = RNG.NormalDistribution(0, 1)
gausian = RNG.CreateGenerator(0, __gausian)

def noiseGen(points, X, Qd, b):
    mhb = -b / 2.0
    Qd = math.sqrt(Qd)  # Deviation of the noise
    hfb = [0] * (points * 2)
    wfb = [0] * (points * 2)
    hfb[0] = 1.0
    wfb[0] = Qd * gausian.ranf()
    for i in range(1, len(wfb)/2):
        # Generate hk coefficients
        hfb[i] = hfb[i-1]/float(i) * (i-1 + mhb)
        # Fille wk with white noise
        wfb[i] = Qd * gausian.ranf()

    hfb = FFT.real_fft(hfb)
    wfb = FFT.real_fft(wfb)

    # Multiply the complex vectors
    # Convolation
    wfb = wfb * hfb
    wfb = FFT.inverse_real_fft(wfb)

    for i in range(0,len(wfb)/2):
        X[i] += wfb[i]


if __name__ == '__main__':
    X = [0] * (2**7)
    noiseGen(2**7, X, 1, -4)
    c = 0
    for i in X:
        print c ,i
        c += 1


> Hi Johan,
>
> > I'm not expecting anybody to look at the whole programs, so I have just
cut
> > out the important part (however, the complete source is included at the
> > bottom of this mail.  The program is creating colored noise)
>
> As a side note, this is a pretty neat function.  Can you give
> me a reference for it?  I would like to know exactly what is
> going on...
>
> > The problem is (most likely) that the C program uses a library called
> > "Numerical Recipes".  In this library there is a function called
realft().
> > I don't know these FFT functions all that well, but there is a distinct
> > difference from the NR (Numerical Recipies) realft() and the
FFT.real_fft()
> > function of Numerical Python.   This is best illustrated by an example:
> > Assume a list/array of 1024 integers.  If you put this array through
> > FFT.real_fft() you get a 513 long array as a result.  The NR realft()
gives
> > a 2048 long array.  Now,  since C cannot deal with complex numbers it
has to
> > use two entries for each number.  Still the array is much larger than
the
> > Numpy version.
> >
> > Anybody know why ?
>
> Well, the FFT module returns an array that contains only the
> positive frequencies (from 0 freq (i.e. the DC value) to the
> Nyquist) from the real-valued FFT.  This is N/2+1 values if N
> is the number of points in the real-valued FFT.
>
> The Numerical Recipes (NR) routine should return N/2 values
> (although you actually get _N_ floats back instead of N/2
> complex numbers -- these are the real and imaginary
> components).  NR also packs the Nyquist and DC components into
> the first two floats (i.e. the first complex number).  This way
> you still get all the information, but in the same number of
> floats as the original array.  If this is confusing, I
> recommend reading the section on how NR packs its FFT arrays.
> The book can be found on-line at:
>
> http://www.ulib.org/webRoot/Books/Numerical_Recipes/bookcpdf.html
>
> > wfb and hfb are equal length. Is this a legal way to convolve using
> > FFT.real_fft() ?
> >
> > wfb = FFT.real_fft(wfb)
> > hfb = FFT.real_fft(hfb)
> >
> > for i in range(0, len(wfb)):
> >        wfb[i] = wfb[i] * hfb[i]
> >
> > wfb = FFT.inverse_real_fft(wfb)
>
> Yes.  But it is not very efficient because of the for loop.  I
> have modified your code to make it array-based (i.e. using the
> wonderful features of Numeric).  Notice that all the for loops
> are gone (or at least hidden somewhere underneath in the C
> implementation...).  I _think_ it does what you want it to do,
> and the only not-quite-so-intuitive thing is the
> umath.multiply.accumulate call, which performs th
> recurrence-like multiplications in the hfb for-loop.
>
> Cheers,
>
> Scott
>
> ---------------
>
> import umath, Numeric, FFT, RandomArray, sys
>
> def noiseGen(points, Qd, b):
>     mhb = -b/2.0
>     Qd = umath.sqrt(Qd)  # Deviation of the noise
>     hfb = Numeric.ones(points, 'd')
>     indices = Numeric.arange(len(hfb)-1, typecode='d')
>     hfb[1:] = (mhb+indices)/(indices+1.0)
>     hfb = umath.multiply.accumulate(hfb)
>     wfb = Qd*RandomArray.standard_normal(points)
>     return FFT.inverse_real_fft(FFT.real_fft(wfb)*FFT.real_fft(hfb))
>
> if __name__ == '__main__':
>     X = noiseGen(2**5, 1, -2)
>     for x in X:  print x
>
> ---------------
>
> --
> Scott M. Ransom              Address:  McGill Univ. Physics Dept.
> Phone:  (514) 398-6492                 3600 University St., Rm 338
> email:  ra...@ph...       Montreal, QC  Canada H3A 2T8
> GPG Fingerprint: 06A9 9553 78BE 16DB 407B  FFCA 9BFA B6FF FFD3 2989
>
> _______________________________________________________________
>
> Have big pipes? SourceForge.net is looking for download mirrors. We supply
> the hardware. You get the recognition. Email Us: ban...@so...
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>

Re: [Numpy-discussion] real_fft

[Scott R. <ra...@ph...>]

Hi Johan,

> I'm not expecting anybody to look at the whole programs, so I have just cut
> out the important part (however, the complete source is included at the
> bottom of this mail.  The program is creating colored noise)

As a side note, this is a pretty neat function.  Can you give
me a reference for it?  I would like to know exactly what is
going on...

> The problem is (most likely) that the C program uses a library called
> "Numerical Recipes".  In this library there is a function called realft().
> I don't know these FFT functions all that well, but there is a distinct
> difference from the NR (Numerical Recipies) realft() and the FFT.real_fft()
> function of Numerical Python.   This is best illustrated by an example:
> Assume a list/array of 1024 integers.  If you put this array through
> FFT.real_fft() you get a 513 long array as a result.  The NR realft() gives
> a 2048 long array.  Now,  since C cannot deal with complex numbers it has to
> use two entries for each number.  Still the array is much larger than the
> Numpy version.
> 
> Anybody know why ?

Well, the FFT module returns an array that contains only the
positive frequencies (from 0 freq (i.e. the DC value) to the
Nyquist) from the real-valued FFT.  This is N/2+1 values if N
is the number of points in the real-valued FFT.

The Numerical Recipes (NR) routine should return N/2 values
(although you actually get _N_ floats back instead of N/2
complex numbers -- these are the real and imaginary
components).  NR also packs the Nyquist and DC components into
the first two floats (i.e. the first complex number).  This way
you still get all the information, but in the same number of
floats as the original array.  If this is confusing, I
recommend reading the section on how NR packs its FFT arrays.
The book can be found on-line at:

http://www.ulib.org/webRoot/Books/Numerical_Recipes/bookcpdf.html

> wfb and hfb are equal length. Is this a legal way to convolve using
> FFT.real_fft() ?
> 
> wfb = FFT.real_fft(wfb)
> hfb = FFT.real_fft(hfb)
> 
> for i in range(0, len(wfb)):
>        wfb[i] = wfb[i] * hfb[i]
> 
> wfb = FFT.inverse_real_fft(wfb)

Yes.  But it is not very efficient because of the for loop.  I
have modified your code to make it array-based (i.e. using the
wonderful features of Numeric).  Notice that all the for loops
are gone (or at least hidden somewhere underneath in the C
implementation...).  I _think_ it does what you want it to do,
and the only not-quite-so-intuitive thing is the
umath.multiply.accumulate call, which performs th
recurrence-like multiplications in the hfb for-loop.

Cheers,

Scott

---------------

import umath, Numeric, FFT, RandomArray, sys

def noiseGen(points, Qd, b):
    mhb = -b/2.0
    Qd = umath.sqrt(Qd)  # Deviation of the noise
    hfb = Numeric.ones(points, 'd')
    indices = Numeric.arange(len(hfb)-1, typecode='d')
    hfb[1:] = (mhb+indices)/(indices+1.0)
    hfb = umath.multiply.accumulate(hfb)
    wfb = Qd*RandomArray.standard_normal(points)
    return FFT.inverse_real_fft(FFT.real_fft(wfb)*FFT.real_fft(hfb))

if __name__ == '__main__':
    X = noiseGen(2**5, 1, -2)
    for x in X:  print x

---------------

-- 
Scott M. Ransom              Address:  McGill Univ. Physics Dept.
Phone:  (514) 398-6492                 3600 University St., Rm 338
email:  ra...@ph...       Montreal, QC  Canada H3A 2T8 
GPG Fingerprint: 06A9 9553 78BE 16DB 407B  FFCA 9BFA B6FF FFD3 2989

[Numpy-discussion] real_fft

[<jo...@oh...>]

Hi there.
I have a problem related to the fft routine in numpy.
I'm trying to port a short C program to python, however it has turned out to
be a little more complicated than initally thought.

I'm not expecting anybody to look at the whole programs, so I have just cut
out the important part (however, the complete source is included at the
bottom of this mail.  The program is creating colored noise)



The problem is (most likely) that the C program uses a library called
"Numerical Recipes".  In this library there is a function called realft().
I don't know these FFT functions all that well, but there is a distinct
difference from the NR (Numerical Recipies) realft() and the FFT.real_fft()
function of Numerical Python.   This is best illustrated by an example:
Assume a list/array of 1024 integers.  If you put this array through
FFT.real_fft() you get a 513 long array as a result.  The NR realft() gives
a 2048 long array.  Now,  since C cannot deal with complex numbers it has to
use two entries for each number.  Still the array is much larger than the
Numpy version.

Anybody know why ?



wfb and hfb are equal length. Is this a legal way to convolve using
FFT.real_fft() ?



wfb = FFT.real_fft(wfb)
hfb = FFT.real_fft(hfb)

for i in range(0, len(wfb)):
       wfb[i] = wfb[i] * hfb[i]

wfb = FFT.inverse_real_fft(wfb)



--
JFØ



====================
PYTHON
====================
import math, FFT, RNG, sys

__gausian   = RNG.NormalDistribution(0, 1)
gausian = RNG.CreateGenerator(0, __gausian)

def noiseGen(points, X, Qd, b):
    mhb = -b / 2.0
    Qd = math.sqrt(Qd)  # Deviation of the noise
    hfb = [0] * points
    wfb = [0] * points
    hfb[0] = 1.0
    wfb[0] = Qd * gausian.ranf()
    for i in range(1, len(wfb)):
        # Generate hk coefficients
        hfb[i] = hfb[i-1]/float(i) * (i-1 + mhb)
        # Fille wk with white noise
        wfb[i] = Qd * gausian.ranf()


    wfb = FFT.real_fft(wfb)
    hfb = FFT.real_fft(hfb)
    print hfb
    # Multiply the complex vectors
    # Convolation
    for i in range(0, len(wfb)):
       wfb[i] = wfb[i] * hfb[i]
    wfb = FFT.inverse_real_fft(wfb)

    for i in range(0,len(wfb)):
        X[i] += wfb[i]


if __name__ == '__main__':
    X = [0] * (2**10)
    noiseGen(2**10, X, 1, -2)
    for i in X:
        print i

==========================
C Program
==========================
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "c/numrec.h"
#include "c/NRutil.h"


void f_beta(int n_pts, float X[], float Q_d, float b, int *idum){
  int i,nn;
  float *hfb,*wfb;
  float mhb,wr,wi;
  nn = n_pts * 2;
  mhb = -b / 2.0;
  Q_d = sqrt(Q_d);

  hfb = vector(1,nn);
  wfb = vector(1,nn);
  hfb[1] = 1.0;
  wfb[1] = Q_d * gasdev(idum);
  for(i=2; i<=n_pts; i++){
    hfb[i] = hfb[i-1]*(mhb+(float)(i-2))/((float)(i-1));
    wfb[i] = Q_d * gasdev(idum);
  }
  for(i=n_pts;i <= nn; i++){
    hfb[i] = 0.0;
    wfb[i] = 0.0;
  }

  realft(hfb, n_pts,1);
  realft(wfb, n_pts,1);

  wfb[1]=wfb[1]*hfb[1];
  wfb[2]=wfb[2]*hfb[2];
  for(i=3; i<= nn; i += 2){
    wr = wfb[i];
    wi = wfb[i+1];
    wfb[i] = wr*hfb[i]-wi*hfb[i+1];
    wfb[i+1] = wr*hfb[i+1]+wi*hfb[i];
  }

  realft(wfb, n_pts, -1);
  for(i=1; i <=n_pts;i++){
    X[i] += wfb[i]/((float)n_pts);
  }
  free_vector(hfb,1,nn);
  free_vector(wfb,1,nn);
}

int main(){
  int length22;
  int i;
  int idum = -210310212;
  float  *X;
  X = vector(1,1024);
  f_beta(1024, X, 1, -2, &idum);
  for(i=1;i<=1024;i++){
    printf("%f\n",X[i]);
  }
  return 0;
}
--
Johan Fredrik Øhman

[Numpy-discussion] tiny bug in MLab.py

[Charles G W. <cg...@al...>]

>>> import Numeric
>>> import MLab
>>> 
>>> a=Numeric.ones(4)
>>> b=Numeric.ones(4)
>>> 
>>> print MLab.corrcoef(a,b)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File "MLab.py", line 230, in corrcoef
    c = cov(x, y)
  File "MLab.py", line 237, in cov
    if y != None: m = concatenate((m,y))
TypeError: function not supported for these types, and can't coerce to
supported types


This is happening because of some comparisons of array objects to
None.  You can't say "a != None" for array objects.  (Maybe you should
be able to?)

Anyhow, the following patch corrects the misbehavior:

--- MLab.py.orig	Fri May 10 16:12:52 2002
+++ MLab.py	Fri May 10 16:36:42 2002
@@ -34,7 +34,7 @@
 	"""eye(N, M=N, k=0, typecode=None) returns a N-by-M matrix where the 
 	k-th diagonal is all ones, and everything else is zeros.
 	"""
-	if M == None: M = N
+	if M is None: M = N
 	if type(M) == type('d'): 
 		typecode = M
 		M = N
@@ -46,7 +46,7 @@
 	"""tri(N, M=N, k=0, typecode=None) returns a N-by-M matrix where all
 	the diagonals starting from lower left corner up to the k-th are all ones.
 	"""
-	if M == None: M = N
+	if M is None: M = N
 	if type(M) == type('d'): 
 		typecode = M
 		M = N
@@ -197,7 +197,7 @@
 	"""trapz(y,x=None,axis=-1) integrates y along the given dimension of
 	the data array using the trapezoidal rule.
 	"""
-	if x == None:
+	if x is None:
 		d = 1.0
 	else:
 		d = diff(x,axis=axis)
@@ -234,7 +234,7 @@
 def cov(m,y=None):
 	m = asarray(m)
 	mu = mean(m)
-	if y != None: m = concatenate((m,y))
+	if y is not None: m = concatenate((m,y))
 	sum_cov = 0.0
 	for v in m:
 		sum_cov = sum_cov+multiply.outer(v,v)

Re: [Numpy-discussion] numarray on Cygwin

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

Jochen Küpper wrote:

>On Tue, 07 May 2002 08:57:00 -0400 Todd Miller wrote:
>
>Todd> Jochen Küpper wrote:
>
>>>I am trying to run numarray on Cygwin.
>>>(Where Numeric seems to work fine.)
>>>
>>>In order to compile it I had to change the use of round to rint, as
>>>Cygwin's C library (newlib) doesn't have round.  The patch is
>>>attached.
>>>
>
>Todd> Trying to call a PyCObject is a symptom of an inconsistent
>Todd> installation.
>
>That was easy.  I am not sure why a 
>,----
>| python setup.py build --force
>| python setup.py install --force
>`----
>didn't do it "right" automatically (it definitely /should/), but
>deleting everyting in site-packages and reinstalling works.
>
>
>The remaining question is how to deal with round:
>a) Use rint instead of round everywhere.
>b) We use rint but make sure we are using the correct rounding mode.
>c) We write our own round function using floor or ceil and use that on
>   Cygwin.
>d) We use rint on Cygwin, round everywhere else.
>
>
>Greetings,
>Jochen
>
I think we should define our own function, based on the definition of 
around in
Numeric.  This is most closely related to choice "c", using floor and 
ceil.  It was an
oversight to use round.  

Todd

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