calculating on matrix indices

[Robert Kern]

Brian Blais wrote:
> Hello,
> 
> In my attempt to learn python, migrating from matlab, I have the following problem. 
> Here is what I want to do, (with the wrong syntax):
> 
> from numpy import *
> 
> t=arange(0,20,.1)
> x=zeros(len(t),'f')
> 
> idx=(t>5)
> tau=5
> x[idx]=exp(-t[idx]/tau)  # <---this line is wrong (gives a TypeError)
> 
> #------------------
> 
> what is the best way to replace the wrong line with something that works: replace all 
> of the values of x at the indices idx with exp(-t/tau) for values of t at indices idx?

The traceback tells you exactly what's wrong:

In [7]: x[idx] = exp(-t[idx]/tau)
---------------------------------------------------------------------------
exceptions.TypeError                                 Traceback (most recent call
last)

/Users/kern/<ipython console>

TypeError: array cannot be safely cast to required type

In [8]: t.dtype
Out[8]: dtype('>f8')

In [9]: x.dtype
Out[9]: dtype('>f4')

In [10]: exp(-t[idx]/tau).dtype
Out[10]: dtype('>f8')

The righthand-side is in double precision, but you are trying to stuff it into a
single-precision array. numpy can't do that conversion without losing precision,
so it raises a TypeError. If you really wanted x to be single precision, then
you would explicitly do the conversion on the righthand-side.

In [11]: x[idx] = exp(-t[idx]/tau).astype(float32)

In [12]: x
Out[12]:
array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.36059493,  0.35345468,  0.34645581,  0.33959553,
...
        0.02024191,  0.0198411 ,  0.01944822,  0.01906312,  0.01868564],
dtype=float32)

However, I suspect that you wanted x to be the same precision as t, which
defaults to double precision just like Python.

In [13]: x = zeros(len(t), float)

In [14]: x[idx] = exp(-t[idx]/tau)

In [15]: x
Out[15]:
array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.36059494,  0.35345468,  0.34645581,  0.33959553,
...
        0.02024191,  0.01984109,  0.01944821,  0.01906311,  0.01868564])

BTW, because using arange with floating point steps is somewhat unreliable
because it's floating point. You may want to use linspace() instead.

In [16]: linspace?
Type:           function
Base Class:     <type 'function'>
String Form:    <function linspace at 0x6e21f0>
Namespace:      Interactive
File:
/Library/Frameworks/Python.framework/Versions/2.4/lib/python2.4/site-packages/numpy-0.9.5.2085-py2.4-macosx-10.4-ppc.egg/numpy/lib/function_base.py
Definition:     linspace(start, stop, num=50, endpoint=True, retstep=False)
Docstring:
    Return evenly spaced numbers.

    Return 'num' evenly spaced samples from 'start' to 'stop'.  If
    'endpoint' is True, the last sample is 'stop'. If 'retstep' is
    True then return the step value used.

-- 
Robert Kern
robert.kern at gmail.com

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter