logical statements (in Python and Numeric)

[Fernando Pérez]

Nathan Given wrote:

> I am graphing piecewise functions, and I have a method of getting the
> proper y values in a vector that works in MATLAB... I was trying to get
> it to work in Python, but it wasn't...
> 
> it has to do with a using an and statement...
> 
> 
> here is my code...
> 
> x = arrayrange(-4.,2.+.001,.001, Float)
> y = zeros(len(x),Float)
> 
> #------ Create Piecewise Function ------#
> y = y + (x < -2)*(((x + 2.)**2) + 4)
> y = y + (x >= -2. and x < 0.) * (x + 6.)
> y = y + (x >=0) * (x ** 2. + 6.)
> 

I agree that the above failure is bloody counterintuitive, to the point where 
I'd call it a bug. It's easy to illustrate with:

In [58]: a=arange(12)

In [59]: a > 5 and a < 8
Out[59]: array([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0])

What you need to use is logical_and():

In [60]: logical_and(a>5,a<8)
Out[60]: array([0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0])

You may consider posting to the Numeric mailing list about this. The 'and' 
operator should either not work at all or work correctly. The above is a 
heinous behavior.

But at least you can use the logical_* functions to get around the problem. 
Your code would then read:

In [49]: x = arrayrange(-4.,2.+.001,.001, Float)

In [50]: y = zeros(len(x),Float)

In [51]: y += (x < -2)*(((x + 2.)**2) + 4)

In [52]: y += logical_and(x>=-2,x<0) * (x+6)

In [53]: y += (x >=0) * (x ** 2. + 6.)


Note the '+=' operator, if I recall correctly it works in-place for Numeric 
arrays, which has memory benefits (significant for large arrays).

For common range creation with a closed interval, I wrote the following 
function which I use constantly. Feel free to use it if you like it. The npts 
optional argument is especially useful, since you often care more about the 
particular size of the list rather than the specific value of delta.

def frange(xini,xfin=None,delta=None,**kw):
    """frange([start,] stop[, step, keywords]) -> list of floats

    Return a Numeric array() containing a progression of floats. Similar to
    arange(), but defaults to a closed interval.

    frange(x0, x1) returns [x0, x0+1, x0+2, ..., x1]; start defaults to 0, and
    the endpoint *is included*. This behavior is different from that of
    range() and arange(). This is deliberate, since frange will probably be
    more useful for generating lists of points for function evaluation, and
    endpoints are often desired in this use. The usual behavior of range() can
    be obtained by setting the keyword 'closed=0', in this case frange()
    basically becomes arange().

    When step is given, it specifies the increment (or decrement). All
    arguments can be floating point numbers.

    frange(x0,x1,d) returns [x0,x0+d,x0+2d,...,xfin] where xfin<=x1.

    frange can also be called with the keyword 'npts'. This sets the number of
    points the list should contain (and overrides the value 'step' might have
    been given). arange() doesn't offer this option.

    Examples:
    >>> frange(3)
    array([ 0.,  1.,  2.,  3.])
    >>> frange(3,closed=0)
    array([ 0.,  1.,  2.])
    >>> frange(1,6,2)
    array([1, 3, 5])
    >>> frange(1,6.5,npts=5)
    array([ 1.   ,  2.375,  3.75 ,  5.125,  6.5  ])
    """
    #defaults
    kw.setdefault('closed',1)
    endpoint = kw['closed'] != 0

    # funny logic to allow the *first* argument to be optional (like range())
    # This was modified with a simpler version from a similar frange() found
    # at http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/66472
    if xfin == None:
        xfin = xini + 0.0
        xini = 0.0

    if delta == None:
        delta = 1.0

    # compute # of points, spacing and return final list
    try:
        npts=kw['npts']
        delta=(xfin-xini)/float(npts-endpoint)
    except KeyError:
        npts=int((xfin-xini)/delta+endpoint)

    return arange(npts)*delta+xini


If you want to be able to plot with gnuplot a single array with g.plot(y), or 
y vs x with g.plot(x,y), ask me for a modified version I have, which also 
fixes the problems with hardcopy(). It's part of ipython 
(http://windom.colorado.edu/~fperez/ipython/) which gives you a python shell 
which feels much like that of IDL or Mathematica. 

Cheers,

f.