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

[Numpy-discussion] rank-0 arrays

[Perry G. <pe...@st...>]

Before we implement what we said we would regarding rank-0
arrays in numarray, there became apparent a couple new issues
that didn't really get considered in the first round of 
discussion (at least I don't recall that they did).

To restate the issue: there was a question about whether
an index to an array that identified a single element only
(i.e., not a slice, nor an incomplete index, e.g. x[3] where
x is two dimensional) should return a Python scalar or a
rank-0 array. Currently Numeric is inconsistent on this point.
One usually gets scalars, but on some occasions, rank-0 arrays
are returned. Good arguments are to be had for either alternative.

The primary advantage of returning rank-0 arrays is that they
reduce the need for conditional code checking to see if a result
is a scalar or an array. At the end of the discussion it was
decided to have numarray return rank-0 arrays in all instances of
single item indexing. Since then, a couple potential snags have
arisen. I've already discussed some of these with Paul Dubois
and Eric Jones. I'd like a little wider input before making a
final (or at least experimental) decision.

If we return rank-0 arrays, what should repr return for rank-0
arrays. My initial impression is that the following is highly
undesirable for a interactive session, but maybe it is just me:

>>> x = arange(10)
>>> x[2]
array(2)

We, of course, could arrange __repr__ to return "2" instead,
in other words print the simple scalar for all cases of rank-0
arrays. This would yield the expected output in the above
example. Nevertheless, isn't it violating the intent of repr?
Are there other examples where Python uses repr in a similar,
misleading manner? But perhaps most feel that returning array(2)
is perfectly acceptable and won't annoy users. I am curious
about what people think about this.

The second issue is an efficiency one. Currently numarray uses
Python objects for arrays. If we return rank-0 arrays for
single item indexing, then some naive uses of larger arrays
as sequences may lead to an enormous number of array objects
to be created. True, there will be equivalent means of doing
the same operation that won't result in massive object creations
(such as specifically converting an array to a list, which would 
be done much faster). Is this a serious problem?

These two issues led us to question whether we should indeed
return rank-0 arrays. We can live with either solution. But
we do want to make the right choice. We also know that both
functionalities must exist, e.g., indexing for scalars and 
indexing for rank-0 arrays and we will provide both. The issue
is what indexing syntax returns. One argument is that it is
not a great burden on programmers to use a method (or other
means) to obtain a rank-0 array always if that is important
for the code they are writing and that we should make the
indexing syntax return what most users (especially less expert
ones) intuitively expect (scalars I presume). But others feel
it is just as important for the syntax that a progammer uses
to be as simple as the interactive user expects (instead of
something like x.getindexasarrayalways(2,4,1) [well, with a much
better, and shorter, name])

Do either of these issues change anyone's opinion? If people
still want rank-0 arrays, what should repr do?

Perry Greenfield
 

[Numpy-discussion] Todd's NumPy Feature Request tracker surprise!

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

I was really surprised today to find that Perry Greenfield and I were 
the only two developers listed on the Numeric Feature Requests tracker 
on Source Forge.  Since we work on numarray,  I removed us.  That leaves 
no one.

If you are a NumPy developer who wants to handle Feature Requests for 
Numeric,  you might want to add yourself back onto the list of assignees.  

Todd

-- 
Todd Miller 			jm...@st...
STSCI / SSG

Re: [Numpy-discussion] Logical AND on arrays

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

Alexander Schmolck wrote:
> Tom Transue <tr...@ni...> writes:
> > useful.  It would be nice to map the first expression to do what the first does.
> 
> Impossible. Python's `and` and `or` are shortcircuting operators with a fixed
> meaning. ``A and B`` only evaluates both ``A`` and ``B`` iff both have the
> boolean value 'True' (e.g. ``0 and 0/1`` is OK because the second part never
> gets evaluated), otherwise the first 'False` value is returned.

true. Another option, however is to use the bitwise and aperator, which
is overloaded:

A & B

This generally will work in the cases you are likely to use it, like:

if (A > 5) & (B < 5):

as you'll be &-ing arrays of zeros and ones.

> there is (luckily) no way to change the meaning of `and` and `or` (just
> as one (luckily) can't change what the `if` statement does).

personally, I don't think it's that lucky. I'd rather lose the slight
benifits of the sometimes confusing ability to short-circuit, and get
full operator overloading. I think "and" and "or" are more like "<" than
"if". However, if I had designed a language it would be a pretty
wretched mess, I'm sure.

-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] Logical AND on arrays

[Alexander S. <a.s...@gm...>]

Tom Transue <tr...@ni...> writes:

> If ar1 and ar2 are arrays, are the following two expressions supposed to give
> the same result?
> 
>   (ar1 and ar2)
> 
>   logical_and(ar1,ar2)
> 
> The first form seems to return the value of the second array, which isn't very
Only if the boolean value of the first array is true (i.e. iff not all its
values are nonzero).

> useful.  It would be nice to map the first expression to do what the first does.

Impossible. Python's `and` and `or` are shortcircuting operators with a fixed
meaning. ``A and B`` only evaluates both ``A`` and ``B`` iff both have the
boolean value 'True' (e.g. ``0 and 0/1`` is OK because the second part never
gets evaluated), otherwise the first 'False` value is returned.

While each class is free to define a truth testing method (`__nonzero__`, or
`__len__` if `__nonzero__` doesn't exist; the convention is that empty
containers and zero numbers are false and pretty much everything else is
true), there is (luckily) no way to change the meaning of `and` and `or` (just
as one (luckily) can't change what the `if` statement does).

If this is still a bit unclear, I think this question has been raised before
in this list, so you should be able to find a thread in the archive.

cheers,

alex


-- 
Alexander Schmolck     Postgraduate Research Student
                       Department of Computer Science
                       University of Exeter
A.S...@gm...     http://www.dcs.ex.ac.uk/people/aschmolc/

[Numpy-discussion] Logical AND on arrays

[Tom T. <tr...@ni...>]

If ar1 and ar2 are arrays, are the following two expressions supposed to give
the same result?

  (ar1 and ar2)

  logical_and(ar1,ar2)

The first form seems to return the value of the second array, which isn't very
useful.  It would be nice to map the first expression to do what the first does.

Tom Transue

[Numpy-discussion] Comparison between Numpy/Numarray and Ox

[<sd...@UC...>]

Dear Numpy list,

Has anyone looked at the comparitive speed of Numpy/Numarray vs. Ox?  For 
those that don't know, Ox is an implicitly typed, object oriented matrix 
programming language that is one of the fastest around, and also allows 
linking to C libraries (http://www.nuff.ox.ac.uk/Users/Doornik/index.html).  I 
use Python a lot, but I was wondering what the computational cost would be to 
switch my matrix-y stuff to Numpy.  Ox isn't open source...

Best wishes
Simon

[Numpy-discussion] [Announce] Numerical Python 22.0

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

Numerical Python release 22.0 is at Sourceforge. A Windows installer and
source zip are also available.

Version 22.0 
    a. Changed multiarraymodule functions to accept keywords where
documentation implies
         it through the use of optional variables.
         Specifically in multiarray: zeros, take, transpose, repeat, 
                                    set_string_function,
cross_correlate.
                      in ufuncobject: reduce and accumulate now take
keyword arguments
                                      for the optional axis argument. 
    b. Added support for unsigned shorts 'w' and unsigned ints 'u' 
         -- Travis Oliphant with help from Darren Hart and F. Oliver
Gathmann.
    Increased max permissible iterations in SVD for supplied lapack. --
Dubois
    Recoded RandomArray.randint to try to see if we can work around bug 
         on some platforms. -- Dubois

Version 21.3 June 8, 2002
Fixed bugs:
    [ #557927 ] fixed matrix slice assignment
    [ #552922 ] added check for correct datatype in .astype() method.
                  Created new API PyArray_ValidType to handle this check
                  here as well as in multiarraymodule.c  
    [ #551808 ] fixed segfault with unicode array (Travis O.)
    [ #559511 ] MLab.std now works for axis != 0 (Travis O.)
    [ #542979 ] sum returns exception tuple

    [ #528328 ] true division operators used to return single precision
on division of
                   integers and longs --- now defaults to double
precision (but only on int and 
                       long division --- still single-precision for
ubyte, short, and byte division.
    [ none   ] arange(start, end, step) slightly different near end
points than 
                   start + arange(0, N)*step where N is the length.
    [ none   ] a = zeros(2,'D');  a[0] = array(0.0+0.6j) would not work.
(Assigning a 
                   rank-0 array did not work for CFLOAT_setitem or
CDOUBLE_setitem.
        [ 530688 ] Python crash when transposing array (Walter Moreira)

Version 21.0 March 13, 2002
Fixed bugs:
    [ #482603 ] Memory leak in MA/Numeric/Python
                Reported by Reggie Dugard. Turned out to be 
                *two* memory leaks in one case in a routine in Numeric, 
                array_objectype. (Dubois)
    [ none    ] if vals was a null-array array([]) putmask and put would
                   crash.  Fixed with check.
    [ #469951 ] n = n1[0] gives array which shares dimension of n1
array. 
                  This causes bugs if shape of n1 is changed (n didn't
used
                  to have it's own dimensions array  (Travis Oliphant)
    [ #514588 ] MLab.cov(x,x) != MLab.cov(x) (Travis Oliphant)
    [ #518702 ] segfault when invalid typecode for asarray (Travis
Oliphant)
    [ #497530 ] MA __getitem__ prevents 0 len arrays (Reggie Duggard)
    [ #508363 ] outerproduct of noncontiguous arrays (Martin Wiechert)
    [ #513010 ] memory leak in comparisons (Byran Nollett)
    [ #512223 ] Character typecode not defined (Jochen Kupper)
    [ #500784 ] MLab.py diff error (anonymous, fixed by Dubois)
    [ #503741 ] accuracy of MLab.std(x) (Katsunori Waragai)
    [ #507568 ] overlapping copy a[2:5] = a[3:6]
                Change uses of memcpy to memmove which allows overlaps.
    [ numpy-Patches-499722 ] size of buffer created from array is bad
(Michel Sanner).
    [ #502186 ] a BUG in RandomArray.normal (introduced by last bug fix
in 20.3)      
               (Katsunori Waragai).
      
Fixed errors for Mac (Jack Jensen).

Make rpm's properly, better Windows installers. (Gerard Vermeulen)
    Added files setup.cfg; setup calculates rpm_install.sh to use
current Python.
    New setup.py, eliminate setup_all.py. Use os.path.join everywhere.
Revision in b6
    added file README.RPM, further improvements.

Implement true division operations for Python 2.2. (Bruce Sherwood)
    Note: true division of all integer types results in an array of
floats, 
    not doubles. This decision is arbitrary and there are arguments
either way,
    so users of this new feature should be aware that the decision may 
    change in the future. 

New functions in Numeric; they work on any sequence a that can be
converted to a
Numeric array. Similar change to average in MA. (Dubois)

    def rank (a):
        "Get the rank of a (the number of dimensions, not a matrix
rank)"

    def shape (a):
        "Get the shape of a"

    def size (a, axis=None):
        "Get the number of elements in a, or along a certain axis."
                            
    def average (a, axis=0, weights=None, returned = 0):
        """average(a, axis=0, weights=None)
           Computes average along indicated axis. 
           If axis is None, average over the entire array.
           Inputs can be integer or floating types; result is type
Float.
       
           If weights are given, result is:
               sum(a*weights)/(sum(weights))
           weights must have a's shape or be the 1-d with length the
size
           of a in the given axis. Integer weights are converted to
Float.

           Not supplying weights is equivalent to supply weights that
are
           all 1.
    
           If returned, return a tuple: the result and the sum of the
weights 
           or count of values. The shape of these two results will be
the same.
    
           raises ZeroDivisionError if appropriate when result is
scalar.
           (The version in MA does not -- it returns masked values).
        """

Re: [Numpy-discussion] tensor ops in dim > 3?

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

"Dr. Dmitry Gokhman" <go...@sp...> writes:

> I have a rank three n-dim tensor A.  For two of the axes I want to perform
> v^t A v (a quadratic form with scalar output, where v is a vector). The
> final output should be a vector.  I also need to compute the derivative of
> this with respect to v.  This involves symmetrizing and matrix-vector
> multiplication (2 sym(A)v using two axes of A only, which gives a vector)
> with the final result being a matrix.

Whenever dealing with somewhat more complex operations of this time, I think
it's best to go back to the basic NumPy functionality rather then figuring
out if there happens to be a function that magically does it. In this case,
assuming that the first axis of A is the one that is not summed over:

   sum( sum(A*v[NewAxis, NewAxis, :], -1) * v[NewAxis, :], -1)

The idea is to align v with one of the dimensions of A, then multiply
elementwise and sum over the common axis. Note that the first (inner)
sum leaves a rank-2 array, so for the second multiplication v gets
extended to rank-2 only.

> PS One more dumb question: I just installed the ScientificPython-2.4.1 rpm
> on my reincarnated Mandrake linux machine running python2.2.  Do I need to
> do something to configure it?  My scripts aren't finding things (e.g.
> indexing.py).

If you took the binary RPM from my site, they might not work correctly
with Mandrake, as they were made for RedHat. The source RPM should
work with all RPM-based Linux distributions. There is nothing that needs
configuring with an RPM.

Also note that Scientific is a package, so the correct way to import the
indexing module is

   import Scientific.indexing

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] tensor ops in dim > 3?

[Travis O. <oli...@ee...>]

> I have a rank three n-dim tensor A.  For two of the axes I want to perform
> v^t A v (a quadratic form with scalar output, where v is a vector). The
> final output should be a vector.  I also need to compute the derivative of
> this with respect to v.  This involves symmetrizing and matrix-vector
> multiplication (2 sym(A)v using two axes of A only, which gives a vector)
> with the final result being a matrix.
>

I'm not exactly sure what you mean by a rank three n-dim tensor (in
Numeric the rank is the number of dimensions in the array).

But, I think you can accomplish what you desire in two different ways:

1) If A has N dimensions and you want to perform the reduction over axis
I'll label a1 and a2 (that is a1 is the axis for the v^t*A  sum while a2
is the axis for the A*v sum).  Then I think this should work (if a2 > a1).

ex1 = [Numeric.NewAxis]*(N-1)
ex1[a1] = slice(None)
ex2 = [Numeric.NewAxis]*N
ex2[a2] = slice(None)

Nar = Numeric.add.reduce
result = Nar(v[ex1]*Nar(A*v[ex2],axis=a2),axis=a1)

# I think you need a recent version of Numeric for the axis keyword
#   to be defined here.  Otherwise, just pass a2 and a1 as arguments
#   without the keyword.


2) Using dot (requires transposing the matrix) as dot only operates over
certain dimensions --- I would not try this.

-Travis Oliphant

[Numpy-discussion] tensor ops in dim > 3?

[Dr. D. G. <go...@sp...>]

I have a rank three n-dim tensor A.  For two of the axes I want to perform
v^t A v (a quadratic form with scalar output, where v is a vector). The
final output should be a vector.  I also need to compute the derivative of
this with respect to v.  This involves symmetrizing and matrix-vector
multiplication (2 sym(A)v using two axes of A only, which gives a vector)
with the final result being a matrix.

Is there something elegant I can do short of extending numpy with fortrash
routines? Perhaps something like SciPy's 3-d tensor ops but for n-dim?

Thanks in advance for any tips, - d

PS One more dumb question: I just installed the ScientificPython-2.4.1 rpm
on my reincarnated Mandrake linux machine running python2.2.  Do I need to
do something to configure it?  My scripts aren't finding things (e.g.
indexing.py).

Re: [Numpy-discussion] numarray: RandomArray2

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

Jochen Küpper wrote:

>On Mon, 19 Aug 2002 09:27:35 -0400 Todd Miller wrote:
>
>Todd> Jochen Küpper wrote:
>
>>>Also note the "strange" docstring of range:
>>>,----
>>>| range(...)
>>>|     range([start,] stop[, step]) -> list of integers
>>>`----
>>>pointing straight toward its behavior.
>>>
>
>Ok, reworked the docstrings of module RandomArray2 a bit in order to
>(hopefully) clarify these issues. I have attached a diff, if it is ok
>I'll check it in.
>
Please do.

Todd

>
>>>So well, maybe someone with insight into RNG's can comment on the,
>>>mirroring issue?
>>>
>
>Anybody?
>
>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.
>
>Put a comment about that in the docstring of random_integers.
>
>Greetings,
>Jochen
>
>
>------------------------------------------------------------------------
>
>? Packages/RandomArray2/build
>? Packages/RandomArray2/Lib/ChangeLog
>Index: Packages/RandomArray2/Lib/RandomArray2.py
>===================================================================
>RCS file: /cvsroot/numpy/numarray/Packages/RandomArray2/Lib/RandomArray2.py,v
>retrieving revision 1.3
>diff -u -r1.3 RandomArray2.py
>--- Packages/RandomArray2/Lib/RandomArray2.py	16 Aug 2002 10:44:21 -0000	1.3
>+++ Packages/RandomArray2/Lib/RandomArray2.py	22 Aug 2002 22:11:16 -0000
>@@ -5,6 +5,11 @@
> import math
> from types import *
> 
>+__doc__ = """Random number array generators.
>+
>+This module provides functions to generate arrays of random numbers.
>+"""
>+
> # Extended RandomArray to provide more distributions:
> # normal, beta, chi square, F, multivariate normal,
> # exponential, binomial, multinomial
>@@ -13,8 +18,8 @@
> ArgumentError = "ArgumentError"
> 
> def seed(x=0,y=0):
>-    """seed(x, y), set the seed using the integers x, y; 
>-    Set a random one from clock if  y == 0
>+    """Set the RNG seed using the integers x, y; 
>+    If |y| == 0 set a random one from clock.
>     """
>     if type (x) != IntType or type (y) != IntType :
>         raise ArgumentError, "seed requires integer arguments."
>@@ -30,14 +35,14 @@
> seed()
> 
> def get_seed():
>-    "Return the current seed pair"
>+    """Return the current seed pair"""
>     return ranlib.get_seeds()
> 
> def _build_random_array(fun, args, shape=[]):
>-# Build an array by applying function fun to
>-# the arguments in args, creating an array with
>-# the specified shape.
>-# Allows an integer shape n as a shorthand for (n,).
>+    """Build an array by applying function |fun| to the arguments in
>+    |args|, creating an array with the specified shape.
>+    Allows an integer shape n as a shorthand for (n,).
>+    """
>     if isinstance(shape, IntType): 
>         shape = [shape]
>     if len(shape) != 0:
>@@ -51,20 +56,28 @@
>         return s[0]
> 
> def random(shape=[]):
>-    "random(n) or random([n, m, ...]) returns array of random numbers"
>+    """random(n) or random([n, m, ...])
>+
>+    Returns array of random numbers in the range 0.0 to 1.0.
>+    """
>     return _build_random_array(ranlib.sample, (), shape)
> 
> def uniform(minimum, maximum, shape=[]):
>-    """uniform(minimum, maximum, shape=[]) returns array of given shape of random reals 
>-    in given range"""
>+    """uniform([minimum,], maximum[, shape])
>+
>+    Return array with shape |shape| of random Floats with all values
>+    > minimum and < maximum.    
>+    """
>     return minimum + (maximum-minimum)*random(shape)
> 
> def randint(minimum, maximum=None, shape=[]):
>-    """randint(min, max, shape=[]) = random integers >=min, < max
>-    If max not given, random integers >= 0, <min"""
>+    """randint([minimum,] maximum[, shape])
>+
>+    Return random integers >= mininimum, < maximum,
>+    if maximum not given, random integers >= 0, < minimum.
>+    """
>     if maximum is None:
>-        maximum = minimum
>-        minimum = 0
>+        maximum, minimum = minimum, 0
>     a = Numeric.floor(uniform(minimum, maximum, shape))
>     if isinstance(a, Numeric.ArrayType):
>         return a.astype(Numeric.Int32)
>@@ -72,79 +85,94 @@
>         return int(a)
>      
> def random_integers(maximum, minimum=1, shape=[]):
>-    """random_integers(max, min=1, shape=[]) = random integers in range min-max inclusive"""
>+    """Return array of random integers in interval [minimum, maximum]
>+
>+    Note that this function has reversed arguments. It is simply a
>+    redirection through randint, and use of that function (randint) is
>+    suggested.
>+    """
>     return randint(minimum, maximum+1, shape) 
>      
> def permutation(n):
>-    "permutation(n) = a permutation of indices range(n)"
>+    """A permutation of indices range(n)"""
>     return Numeric.argsort(random(n))
> 
> def standard_normal(shape=[]):
>-    """standard_normal(n) or standard_normal([n, m, ...]) returns array of
>-           random numbers normally distributed with mean 0 and standard
>-           deviation 1"""
>+    """standard_normal(n) or standard_normal([n, m, ...])
>+
>+    Returns array of random numbers normally distributed with mean 0
>+    and standard deviation 1.
>+    """
>     return _build_random_array(ranlib.standard_normal, (), shape)
> 
> def normal(mean, std, shape=[]):
>-        """normal(mean, std, n) or normal(mean, std, [n, m, ...]) returns
>-           array of random numbers randomly distributed with specified mean and
>-           standard deviation"""
>-        s = standard_normal(shape)
>-        return s * std + mean
>+    """normal(mean, std, n) or normal(mean, std, [n, m, ...])
>+
>+    Returns array of random numbers randomly distributed with
>+    specified mean and standard deviation
>+    """
>+    s = standard_normal(shape)
>+    return s * std + mean
> 
> def multivariate_normal(mean, cov, shape=[]):
>-       """multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...])
>-          returns an array containing multivariate normally distributed random numbers
>-          with specified mean and covariance.
>-
>-          mean must be a 1 dimensional array. cov must be a square two dimensional
>-          array with the same number of rows and columns as mean has elements.
>-
>-          The first form returns a single 1-D array containing a multivariate
>-          normal.
>-
>-          The second form returns an array of shape (m, n, ..., cov.getshape()[0]).
>-          In this case, output[i,j,...,:] is a 1-D array containing a multivariate
>-          normal."""
>-       # Check preconditions on arguments
>-       mean = Numeric.array(mean)
>-       cov = Numeric.array(cov)
>-       if len(mean.getshape()) != 1:
>-              raise ArgumentError, "mean must be 1 dimensional."
>-       if (len(cov.getshape()) != 2) or (cov.getshape()[0] != cov.getshape()[1]):
>-              raise ArgumentError, "cov must be 2 dimensional and square."
>-       if mean.getshape()[0] != cov.getshape()[0]:
>-              raise ArgumentError, "mean and cov must have same length."
>-       # Compute shape of output
>-       if isinstance(shape, IntType): shape = [shape]
>-       final_shape = list(shape[:])
>-       final_shape.append(mean.getshape()[0])
>-       # Create a matrix of independent standard normally distributed random
>-       # numbers. The matrix has rows with the same length as mean and as
>-       # many rows are necessary to form a matrix of shape final_shape.
>-       x = ranlib.standard_normal(Numeric.multiply.reduce(final_shape))
>-       x.setshape(Numeric.multiply.reduce(final_shape[0:len(final_shape)-1]),
>-                  mean.getshape()[0])
>-       # Transform matrix of standard normals into matrix where each row
>-       # contains multivariate normals with the desired covariance.
>-       # Compute A such that matrixmultiply(transpose(A),A) == cov.
>-       # Then the matrix products of the rows of x and A has the desired
>-       # covariance. Note that sqrt(s)*v where (u,s,v) is the singular value
>-       # decomposition of cov is such an A.
>-       (u,s,v) = LinearAlgebra2.singular_value_decomposition(cov)
>-       x = Numeric.matrixmultiply(x*Numeric.sqrt(s),v)
>-       # The rows of x now have the correct covariance but mean 0. Add
>-       # mean to each row. Then each row will have mean mean.
>-       Numeric.add(mean,x,x)
>-       x.setshape(final_shape)
>-       return x
>+    """multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...])
>+
>+    Returns an array containing multivariate normally distributed
>+    random numbers with specified mean and covariance.
>+
>+    |mean| must be a one-dimensional array. |cov| must be a square
>+    two-dimensional array with the same number of rows and columns as
>+    |mean| has elements.
>+
>+    The first form returns a single 1-D array containing a
>+    multivariate normal.
>+
>+    The second form returns an array of shape (m, n, ...,
>+    cov.getshape()[0]). In this case, output[i,j,...,:] is a 1-D array
>+    containing a multivariate normal.
>+    """
>+    # Check preconditions on arguments
>+    mean = Numeric.array(mean)
>+    cov = Numeric.array(cov)
>+    if len(mean.getshape()) != 1:
>+        raise ArgumentError, "mean must be 1 dimensional."
>+    if (len(cov.getshape()) != 2) or (cov.getshape()[0] != cov.getshape()[1]):
>+        raise ArgumentError, "cov must be 2 dimensional and square."
>+    if mean.getshape()[0] != cov.getshape()[0]:
>+        raise ArgumentError, "mean and cov must have same length."
>+    # Compute shape of output
>+    if isinstance(shape, IntType): shape = [shape]
>+    final_shape = list(shape[:])
>+    final_shape.append(mean.getshape()[0])
>+    # Create a matrix of independent standard normally distributed random
>+    # numbers. The matrix has rows with the same length as mean and as
>+    # many rows are necessary to form a matrix of shape final_shape.
>+    x = ranlib.standard_normal(Numeric.multiply.reduce(final_shape))
>+    x.setshape(Numeric.multiply.reduce(final_shape[0:len(final_shape)-1]),
>+               mean.getshape()[0])
>+    # Transform matrix of standard normals into matrix where each row
>+    # contains multivariate normals with the desired covariance.
>+    # Compute A such that matrixmultiply(transpose(A),A) == cov.
>+    # Then the matrix products of the rows of x and A has the desired
>+    # covariance. Note that sqrt(s)*v where (u,s,v) is the singular value
>+    # decomposition of cov is such an A.
>+    (u,s,v) = LinearAlgebra2.singular_value_decomposition(cov)
>+    x = Numeric.matrixmultiply(x*Numeric.sqrt(s),v)
>+    # The rows of x now have the correct covariance but mean 0. Add
>+    # mean to each row. Then each row will have mean mean.
>+    Numeric.add(mean,x,x)
>+    x.setshape(final_shape)
>+    return x
> 
> def exponential(mean, shape=[]):
>-    """exponential(mean, n) or exponential(mean, [n, m, ...]) returns array
>-      of random numbers exponentially distributed with specified mean"""
>-   # If U is a random number uniformly distributed on [0,1], then
>-   #      -ln(U) is exponentially distributed with mean 1, and so
>-   #      -ln(U)*M is exponentially distributed with mean M.
>+    """exponential(mean, n) or exponential(mean, [n, m, ...])
>+
>+    Returns array of random numbers exponentially distributed with
>+    specified mean
>+    """
>+    # If U is a random number uniformly distributed on [0,1], then
>+    #      -ln(U) is exponentially distributed with mean 1, and so
>+    #      -ln(U)*M is exponentially distributed with mean M.
>     x = random(shape)
>     Numeric.log(x, x)
>     Numeric.subtract(0.0, x, x)
>@@ -160,52 +188,79 @@
>     return _build_random_array(ranlib.gamma, (a, r), shape)
> 
> def F(dfn, dfd, shape=[]):
>-    """F(dfn, dfd) or F(dfn, dfd, [n, m, ...]) returns array of F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator."""
>+    """F(dfn, dfd) or F(dfn, dfd, [n, m, ...])
>+
>+    Returns array of F distributed random numbers with dfn degrees of
>+    freedom in the numerator and dfd degrees of freedom in the
>+    denominator.
>+    """
>     return ( chi_square(dfn, shape) / dfn) / ( chi_square(dfd, shape) / dfd)
> 
> def noncentral_F(dfn, dfd, nconc, shape=[]):
>-    """noncentral_F(dfn, dfd, nonc) or noncentral_F(dfn, dfd, nonc, [n, m, ...]) returns array of noncentral F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator, and noncentrality parameter nconc."""
>+    """noncentral_F(dfn, dfd, nonc) or noncentral_F(dfn, dfd, nonc, [n, m, ...])
>+
>+    Returns array of noncentral F distributed random numbers with dfn
>+    degrees of freedom in the numerator and dfd degrees of freedom in
>+    the denominator, and noncentrality parameter nconc.
>+    """
>     return ( noncentral_chi_square(dfn, nconc, shape) / dfn ) / ( chi_square(dfd, shape) / dfd )
> 
> def chi_square(df, shape=[]):
>-    """chi_square(df) or chi_square(df, [n, m, ...]) returns array of chi squared distributed random numbers with df degrees of freedom."""
>+    """chi_square(df) or chi_square(df, [n, m, ...])
>+
>+    Returns array of chi squared distributed random numbers with df
>+    degrees of freedom.
>+    """
>     return _build_random_array(ranlib.chisquare, (df,), shape)
> 
> def noncentral_chi_square(df, nconc, shape=[]):
>-    """noncentral_chi_square(df, nconc) or chi_square(df, nconc, [n, m, ...]) returns array of noncentral chi squared distributed random numbers with df degrees of freedom and noncentrality parameter."""
>+    """noncentral_chi_square(df, nconc) or chi_square(df, nconc, [n, m, ...])
>+
>+    Returns array of noncentral chi squared distributed random numbers
>+    with df degrees of freedom and noncentrality parameter.
>+    """
>     return _build_random_array(ranlib.noncentral_chisquare, (df, nconc), shape)
> 
> def binomial(trials, p, shape=[]):
>-    """binomial(trials, p) or binomial(trials, p, [n, m, ...]) returns array of binomially distributed random integers.
>+    """binomial(trials, p) or binomial(trials, p, [n, m, ...])
>+
>+    Returns array of binomially distributed random integers.
> 
>-           trials is the number of trials in the binomial distribution.
>-           p is the probability of an event in each trial of the binomial distribution."""
>+    |trials| is the number of trials in the binomial distribution.
>+    |p| is the probability of an event in each trial of the binomial
>+    distribution.
>+    """
>     return _build_random_array(ranlib.binomial, (trials, p), shape)
> 
> def negative_binomial(trials, p, shape=[]):
>-    """negative_binomial(trials, p) or negative_binomial(trials, p, [n, m, ...]) returns
>-           array of negative binomially distributed random integers.
>+    """negative_binomial(trials, p) or negative_binomial(trials, p, [n, m, ...])
>+
>+    Returns array of negative binomially distributed random integers.
> 
>-           trials is the number of trials in the negative binomial distribution.
>-           p is the probability of an event in each trial of the negative binomial distribution."""
>+    |trials| is the number of trials in the negative binomial
>+    distribution. |p| is the probability of an event in each trial of
>+    the negative binomial distribution.
>+    """
>     return _build_random_array(ranlib.negative_binomial, (trials, p), shape)
> 
> def multinomial(trials, probs, shape=[]):
>-    """multinomial(trials, probs) or multinomial(trials, probs, [n, m, ...]) returns
>-           array of multinomial distributed integer vectors.
>+    """multinomial(trials, probs) or multinomial(trials, probs, [n, m, ...])
>+
>+    Returns array of multinomial distributed integer vectors.
>+
>+    |trials| is the number of trials in each multinomial distribution. 
>+    |probs| is a one dimensional array. There are len(prob)+1 events. 
>+    prob[i] is the probability of the i-th event, 0<=i<len(prob). The
>+    probability of event len(prob) is 1.-Numeric.sum(prob).
> 
>-           trials is the number of trials in each multinomial distribution.
>-           probs is a one dimensional array. There are len(prob)+1 events. 
>-           prob[i] is the probability of the i-th event, 0<=i<len(prob).
>-           The probability of event len(prob) is 1.-Numeric.sum(prob).
>-
>-       The first form returns a single 1-D array containing one multinomially
>-           distributed vector.
>-
>-           The second form returns an array of shape (m, n, ..., len(probs)).
>-           In this case, output[i,j,...,:] is a 1-D array containing a multinomially
>-           distributed integer 1-D array."""
>-        # Check preconditions on arguments
>+    The first form returns a single 1-D array containing one
>+    multinomially distributed vector.
>+
>+    The second form returns an array of shape (m, n, ..., len(probs)). 
>+    In this case, output[i,j,...,:] is a 1-D array containing a
>+    multinomially distributed integer 1-D array.
>+    """
>+    # Check preconditions on arguments
>     probs = Numeric.array(probs)
>     if len(probs.getshape()) != 1:
>         raise ArgumentError, "probs must be 1 dimensional."
>@@ -215,14 +270,18 @@
>     final_shape.append(probs.getshape()[0]+1)
>     x = ranlib.multinomial(trials, probs.astype(Numeric.Float32),
>                            Numeric.multiply.reduce(shape))
>-        # Change its shape to the desire one
>+    # Change its shape to the desire one
>     x.setshape(final_shape)
>     return x
> 
> def poisson(mean, shape=[]):
>-    """poisson(mean) or poisson(mean, [n, m, ...]) returns array of poisson
>-           distributed random integers with specifed mean."""
>+    """poisson(mean) or poisson(mean, [n, m, ...])
>+
>+    Returns array of poisson distributed random integers with specifed
>+    mean.
>+    """
>     return _build_random_array(ranlib.poisson, (mean,), shape)
>+
> 
> def test():
>     import test as _test
>Index: Packages/RandomArray2/Lib/__init__.py
>===================================================================
>RCS file: /cvsroot/numpy/numarray/Packages/RandomArray2/Lib/__init__.py,v
>retrieving revision 1.1
>diff -u -r1.1 __init__.py
>--- Packages/RandomArray2/Lib/__init__.py	21 Jun 2002 18:25:29 -0000	1.1
>+++ Packages/RandomArray2/Lib/__init__.py	22 Aug 2002 22:11:16 -0000
>@@ -2,3 +2,6 @@
> 
> from RandomArray2 import *
> 
>+__doc__ = RandomArray2.__doc__ + """
>+See RandomArray2.RandomArray2 for more information.
>+"""
>


-- 
Todd Miller 			jm...@st...
STSCI / SSG

Re: [Numpy-discussion] numarray: RandomArray2

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

On Mon, 19 Aug 2002 09:27:35 -0400 Todd Miller wrote:

Todd> Jochen K=FCpper wrote:

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

Ok, reworked the docstrings of module RandomArray2 a bit in order to
(hopefully) clarify these issues. I have attached a diff, if it is ok
I'll check it in.

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

Anybody?

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.

Put a comment about that in the docstring of random_integers.

Greetings,
Jochen
--=20
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] numarray, Convolve, and lineshapes

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

Jochen Küpper wrote:

>-----BEGIN PGP SIGNED MESSAGE-----
>Hash: SHA1
>
>On Thu, 22 Aug 2002 14:25:37 -0400 Todd Miller wrote:
>
>Todd> Your contribution to the Convolve package is a welcome addition
>Todd> provided that you are willing to tighten up your copyrights a
>Todd> little.
>
>That's fine, no problem.
>
Excellent.

>
>
>You can assume that I am willing to allow a 'modified BSD'-style
>license [1] on everything that I'll contribute to numarray. I might
>
Good.

>
>still send it to the ml as GPL'ed, because that's what I put all my
>code under initially. (And if it doesn't go into numarray, I'd
>
That's OK.

>
>probably like to keep it that way.)
>
>Todd> open source under a BSD-like license.
>
My mistake here.   numarray is already using a modified-BSD-like 
license.    

>
>
>I don't like "normal" BSD so much because I wish to release most of my
>code GPL'ed, unless somebody asks. But GPL and "normal" BSD don't get
>along well. (See [1])
>
To keep open the possibility that numarray might someday be included in 
Python, we don't want any part of numarray to be GPL'ed.   This is a 
little overkill, because Packages probably wouldn't be on the immediate 
list of things to submit anyway.   Still,  it keeps life simple and us 
out of trouble.

>
>So I'll clean it up a little tonight and check it in, if that's ok
>with you.
>
Sounds great!

>
>Greetings,
>Jochen
>
>Footnotes: 
>[1]  http://www.gnu.org/philosophy/license-list.html#ModifiedBSD
>- -- 
>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
>-----BEGIN PGP SIGNATURE-----
>Version: GnuPG v1.0.6-cygwin-fcn-1 (Cygwin)
>Comment: Processed by Mailcrypt and GnuPG <http://www.gnupg.org/>
>
>iEYEARECAAYFAj1lMK0ACgkQiJ/aUUS8zY7fmgCfd+3XlZSkdGT/tfOpEZEh3bmu
>jRUAni3NHLtJWsIumxpkO6anLUyonyPF
>=EUBk
>-----END PGP SIGNATURE-----
>
>
>
>-------------------------------------------------------
>This sf.net email is sponsored by: OSDN - Tired of that same old
>cell phone?  Get a new here for FREE!
>https://www.inphonic.com/r.asp?r=sourceforge1&refcode1=vs3390
>_______________________________________________
>Numpy-discussion mailing list
>Num...@li...
>https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>

Regards,
Todd

-- 
Todd Miller 			jm...@st...
STSCI / SSG

Re: [Numpy-discussion] numarray, Convolve, and lineshapes

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

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On Thu, 22 Aug 2002 14:25:37 -0400 Todd Miller wrote:

Todd> Your contribution to the Convolve package is a welcome addition
Todd> provided that you are willing to tighten up your copyrights a
Todd> little.

That's fine, no problem.

You can assume that I am willing to allow a 'modified BSD'-style
license [1] on everything that I'll contribute to numarray. I might
still send it to the ml as GPL'ed, because that's what I put all my
code under initially. (And if it doesn't go into numarray, I'd
probably like to keep it that way.)

Todd> open source under a BSD-like license.

I don't like "normal" BSD so much because I wish to release most of my
code GPL'ed, unless somebody asks. But GPL and "normal" BSD don't get
along well. (See [1])

So I'll clean it up a little tonight and check it in, if that's ok
with you.

Greetings,
Jochen

Footnotes: 
[1]  http://www.gnu.org/philosophy/license-list.html#ModifiedBSD
- -- 
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
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.0.6-cygwin-fcn-1 (Cygwin)
Comment: Processed by Mailcrypt and GnuPG <http://www.gnupg.org/>

iEYEARECAAYFAj1lMK0ACgkQiJ/aUUS8zY7fmgCfd+3XlZSkdGT/tfOpEZEh3bmu
jRUAni3NHLtJWsIumxpkO6anLUyonyPF
=EUBk
-----END PGP SIGNATURE-----

Re: [Numpy-discussion] numarray, Convolve, and lineshapes

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

Hi Jochen,

Your contribution to the Convolve package is a welcome addition provided 
that you are willing to tighten up your copyrights a little. 
 Doc/unittest.py has an example of an acceptable copyright, where 
"Python" should technically be replaced by "numarray",  but where the 
intent remains the same:  open source under a BSD-like license.

Regards,
Todd

Jochen Küpper wrote:

>Looking at Convolve last night I came up with the idea that we often
>need to convolute datasets with common profiles, such as Gauss- or
>Voigt-lineshapes, for example.
>
>So tonight I took my old module I had for that purpose apart and put
>it as submodule into Convolve. Well, currently it really only supports
>Gauss, Lorentz, and Voigt, but the general "framework" to add more
>profiles is there. Adding a new one would consist of defining a new
>derived class (overall 5 lines of code) and providing a function that
>calculates the actual line profile.
>
>I chose the implementation using functor classes so profiles can be
>reused.
>
>I am aware that there are some issues with the functions in
>Lineshape.py (error checking, not automatically converting sequences
>to arrays, ...). If this is a way to go I would move the functions to
>C anyway. Nevertheless a sample function in Python will be provided,
>so new profiles can easily added without caring about C.
>
>I would like to get some feedback whether you think this should be
>included in numarray. I believe it fits very well into Convolve.
>
>
>Example use:
>
>import numarray as num
>import random as ran
>import Convolve
>import Convolve.Lineshape as ls
>
>resolution = 0.1
>fwhm = 0.5
>
># get some (random) data, i.e. a stick-spectrum
>x = num.arange(-50, 50+resolution, resolution)
>y = num.zeros(x.shape, num.Float)
>for i in range(10):
>    y[ran.randint(0, len(y)-1)] = ran.random()
>
># create lineshape objects
>g = ls.GaussProfile(num.arange(-10*fwhm, 10*fwhm+resolution, resolution), 1.0)
>v = ls.VoigtProfile(num.arange(-10*fwhm, 10*fwhm+resolution, resolution), (0.6, 0.6))
>
># convolute data with profile
>res_g = Convolve.convolve(y, g(), Convolve.SAME)
>res_v = Convolve.convolve(y, v(), Convolve.SAME)
>for i in zip(x, y, res_g, res_v):
>    print "%12.6f %12.6f %12.6f %12.6f" % i
>
>
>I attach an archive of the whole Convolve/ dir here in my local copy.
>
>Greetings,
>Jochen
>


-- 
Todd Miller 			jm...@st...
STSCI / SSG

[Numpy-discussion] numarray, Convolve, and lineshapes

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

Looking at Convolve last night I came up with the idea that we often
need to convolute datasets with common profiles, such as Gauss- or
Voigt-lineshapes, for example.

So tonight I took my old module I had for that purpose apart and put
it as submodule into Convolve. Well, currently it really only supports
Gauss, Lorentz, and Voigt, but the general "framework" to add more
profiles is there. Adding a new one would consist of defining a new
derived class (overall 5 lines of code) and providing a function that
calculates the actual line profile.

I chose the implementation using functor classes so profiles can be
reused.

I am aware that there are some issues with the functions in
Lineshape.py (error checking, not automatically converting sequences
to arrays, ...). If this is a way to go I would move the functions to
C anyway. Nevertheless a sample function in Python will be provided,
so new profiles can easily added without caring about C.

I would like to get some feedback whether you think this should be
included in numarray. I believe it fits very well into Convolve.


Example use:

import numarray as num
import random as ran
import Convolve
import Convolve.Lineshape as ls

resolution =3D 0.1
fwhm =3D 0.5

# get some (random) data, i.e. a stick-spectrum
x =3D num.arange(-50, 50+resolution, resolution)
y =3D num.zeros(x.shape, num.Float)
for i in range(10):
    y[ran.randint(0, len(y)-1)] =3D ran.random()

# create lineshape objects
g =3D ls.GaussProfile(num.arange(-10*fwhm, 10*fwhm+resolution, resolution),=
 1.0)
v =3D ls.VoigtProfile(num.arange(-10*fwhm, 10*fwhm+resolution, resolution),=
 (0.6, 0.6))

# convolute data with profile
res_g =3D Convolve.convolve(y, g(), Convolve.SAME)
res_v =3D Convolve.convolve(y, v(), Convolve.SAME)
for i in zip(x, y, res_g, res_v):
    print "%12.6f %12.6f %12.6f %12.6f" % i


I attach an archive of the whole Convolve/ dir here in my local copy.

Greetings,
Jochen
--=20
Einigkeit und Recht und Freiheit                http://www.Jochen-Kuepper.de
    Libert=E9, =C9galit=E9, Fraternit=E9                GnuPG key: 44BCCD8E
        Sex, drugs and rock-n-roll

Re: [Numpy-discussion] Re: numarray -- correlate

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

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On Wed, 21 Aug 2002 08:49:48 -0400 Todd Miller wrote:

Todd> Jochen K=FCpper wrote:

>> Is there any good reason to not rename Convolve.cross_correlate to
>> correlate?  That fits the other names much better.
>>=20
Todd> I agree.

>>=20
>> If needed (aka backwards-compatibility to numpy) we could even provide
>> an alias cross_correlate -- probably documented to be deprecated...=20

Do we need this?

>> But since you have to load a separate module anyway, I think we can
>> rename the function just as well.
>>=20
Todd> Please do this.

Done.

Greetings,
Jochen
- --=20
Einigkeit und Recht und Freiheit                http://www.Jochen-Kuepper.de
    Libert=E9, =C9galit=E9, Fraternit=E9                GnuPG key: 44BCCD8E
        Sex, drugs and rock-n-roll
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.0.6 (GNU/Linux)
Comment: Processed by Mailcrypt and GnuPG <http://www.gnupg.org/>

iD8DBQE9Y58ciJ/aUUS8zY4RAqeBAKCeAjeXmvl01zgU5trFmh6hteeF1wCgvYYY
cbmMNmYlVQF5rqPBczBYsq4=3D
=3D6eTq
-----END PGP SIGNATURE-----

[Numpy-discussion] Re: numarray -- correlate

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

Jochen Küpper wrote:

>Is there any good reason to not rename Convolve.cross_correlate to
>correlate?  That fits the other names much better.
>
I agree.

>
>If needed (aka backwards-compatibility to numpy) we could even provide
>an alias cross_correlate -- probably documented to be deprecated... But
>since you have to load a separate module anyway, I think we can rename
>the function just as well.
>
Please do this.

>
>Greetings,
>Jochen
>

Todd

-- 
Todd Miller 			jm...@st...
STSCI / SSG

[Numpy-discussion] numarray -- correlate

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

Is there any good reason to not rename Convolve.cross_correlate to
correlate?  That fits the other names much better.

If needed (aka backwards-compatibility to numpy) we could even provide
an alias cross_correlate -- probably documented to be deprecated... But
since you have to load a separate module anyway, I think we can rename
the function just as well.

Greetings,
Jochen
--=20
Einigkeit und Recht und Freiheit                http://www.Jochen-Kuepper.de
    Libert=E9, =C9galit=E9, Fraternit=E9                GnuPG key: 44BCCD8E
        Sex, drugs and rock-n-roll

[Numpy-discussion] RE: Cygwin floating point control

[Norman V. <nh...@ca...>]

Jochen K=FCpper writes:
>
>The following program does not link on Cygwin
>,----
>| #include <ieeefp.h>
>|
>| int main()
>| {
>|     return(fgetsticky());
>| }
>`----
>due to fpgetsticky beeing unresolved.
>
>fpgetsticky is declared in /usr/include/ieeefp.h, but I cannot find
>any library it is defined in.  A search on the Cywin site and at
>google didn't turn up anything for fpgetsticky on Cygwin. (?)
>
>Is it implemented?  If so, where?  If not, is there any substitute I
>should be using?

I use this 'home spun' header that I 'hacked' from the BSD sources

Not really sure if it is 100% right or if the license allows this to be
used in Cygwin 'proper' but ... it seems to work for me

HTH

Norman

[Numpy-discussion] Cygwin floating point control

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

The following program does not link on Cygwin
,----
| #include <ieeefp.h>
| 
| int main() 
| {
|     return(fgetsticky());
| }
`----
due to fpgetsticky beeing unresolved.

fpgetsticky is declared in /usr/include/ieeefp.h, but I cannot find
any library it is defined in.  A search on the Cywin site and at
google didn't turn up anything for fpgetsticky on Cygwin. (?)

Is it implemented?  If so, where?  If not, is there any substitute I
should be using?

Greetings,
Jochen

PS: It would be nice if you could cc me.
-- 
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

[Numpy-discussion] numarray on Cygwin

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

Hi,

I could succesfully compile numarray (see numpy.org) on Cygwin
(updated last week), using the standard Cygwin python.  However many
tests fail, due to "missing" overflow or "division by zero"
exceptions.

I took a quick look at the test-failures of numarray on numpy. The
first problem is that we don't get a "Overflow error" warning when it
is expected:
,----[Cygwin]
| >>> array([1, 8, 100, 100], type=Int8) * array([1, 8, 100, -100], type=Int8)
| array([   1,   64,  127, -128], type=Int8)
`----

Whereas on Linux we get 
,----[Linux]
| >>> array([1, 8, 100, 100], type=Int8) * array([1, 8, 100, -100], type=Int8)
| Warning: Encountered overflow(s)  in multiply
| array([   1,   64,  127, -128], type=Int8)
`----

Is this related to a configuration option of python (--with-sigfpe ?)
or a "feature" of the Cygwin libm/libc?

Any ideas?


Btw, I tried to rebuild python (cvs maint22-release patched for
_socket) but got the following errors with the math module (similar
for cmath):
,----
| gcc -shared -Wl,--enable-auto-image-base build/temp.cygwin-1.3.12-i686-2.2/mathmodule.o -L/usr/local/lib -L. -lm -lpython2.2 -o build/lib.cygwin-1.3.12-i686-2.2/math.dll
| build/temp.cygwin-1.3.12-i686-2.2/mathmodule.o: In function `math_1':
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:57: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:57: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:57: undefined reference to `PyFPE_jbuf'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:57: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:59: undefined reference to `PyFPE_dummy'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:59: undefined reference to `PyFPE_counter'
| build/temp.cygwin-1.3.12-i686-2.2/mathmodule.o: In function `math_2':
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:74: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:74: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:74: undefined reference to `PyFPE_jbuf'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:74: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:76: undefined reference to `PyFPE_dummy'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:76: undefined reference to `PyFPE_counter'
| build/temp.cygwin-1.3.12-i686-2.2/mathmodule.o: In function `math_ldexp':
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:173: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:173: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:173: undefined reference to `PyFPE_jbuf'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:173: undefined reference to `PyFPE_counter'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:175: undefined reference to `PyFPE_dummy'
| /home/software/programming/compiler/python/dist/src/Modules/mathmodule.c:175: undefined reference to `PyFPE_counter'
| collect2: ld returned 1 exit status
| WARNING: building of extension "math" failed: command 'gcc' failed with exit status 1
`----

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

[Numpy-discussion] Re: numarray

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

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

>> The following patch is required for Cygwin:
>> 
>> Index: Src/_chararraymodule.c
>> 
Todd> I'll apply this.

Applied.

>> ,----
>> | *****************************************************************
>> | Failure in example:
>> | array([1, 8, 100, 100], type=Int8) * array([1, 8, 100, -100], type=Int8)
>> | from line #1671 of numtest
>> | Expected:
>> | Warning: Encountered overflow(s)  in multiply
>> | array([   1,   64,  127, -128], type=Int8)
>> | Got: array([   1,   64,  127, -128], type=Int8)
>> | *****************************************************************
>> `----

Todd> This overflow exception issue came up earlier this year with
Todd> IRIX.  My guess is that the function int_overflow_error() in
Todd> codegenerator.py is not working correctly.  The suggested
Todd> temporary patch at the time was:

[patch against codegenerator.py]

Todd> If you try it out and it solves the problem for Cygwin, we
Todd> should probably change it universally.

No, it doesn't solve the problem.  I hope to find some time later this
week to look at it.

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

[Numpy-discussion] Re: numarray: RandomArray2

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

Jochen Küpper wrote:

>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"  :))
>
I read it more like:  someone wanted to make one particular user happy 
by giving them the interface they asked for.   Since writing 
random_integers was easy, they did it even though it amounted to a small 
duplication of effort and interface functionality.  But, I'm just 
speculating.

>
>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)
>`----
>??? 
>
Perhaps my intuition about "range" was incorrect, or we're taking the 
analogy too far.  It seems to me that randint is returning random 
numbers between two limits, and the order of the limits is irrelevant. 
 randint(2,-4)  is identical to randint(-4,2). As you have pointed out, 
this is distinctly different than range.   At this point, I'm still 
unconvinced that the world would be made a better place either by 
raising an exception or by changing the semantics to return an "empty 
array".  That small change seems just as likely to upset some user as to 
make another happy.  

>Your-even-more-confused-ly's,
>Jochen
>

Todd
-- 
Todd Miller 			jm...@st...
STSCI / SSG

[Numpy-discussion] Pyrex now interfaces with NumPy (old) array internals

[Frank H. <fra...@cs...>]

Perhaps this is best left to Greg Ward, but I just thought that I'd
point out to the list that Pyrex
<http://www.cosc.canterbury.ac.nz/~greg/python/Pyrex/>
in the latest version 0.4.1 successfully interfaces with the (original)
Numeric array data types. There are examples in both the Demo
subdirectory and in the website's FAQ.

Another alternative for getting performance, interfacing to C libraries,
and even calling Fortran (from Pyrex's Python/C hybrid code).

Frank (Choice is Good ;-) Horowitz