[SciPy-user] numarray complex comparisons

[Tim Hochberg]

[Eric]
[Some stuff elided
> I'll state the case on last time:
>
> One of the beauties of array math is that algorithms are clear and
> concise without a lot of control gunk to obscure what is going on.  If
> exceptions are thrown, every comparison in a "generic algorithm" in
> SciPy will have to be wrapped with a try/except clause or have a lot of
> pre-test code to check for complex values.  The same is true for minimum
> and maximum calculations.  As a quick example:
>
> def clip(a,low,hi):
>     a = asarray(a) # allows lists or scalars as input
>     res = where(a < low,low,a)
>     res = where(res > hi,hi,res)
>     return res
>
> def clip2(a,low,hi):
>     a = asarray(a)
>     if a.typecode() in ['F','D']:
>         reals = a.real
>     else:
>         reals = a
>     if low is complex:
>         rlow = low.real
>     else:
>         rlow = low
>     if hi is complex:
>         rhi = hi.real
>     else:
>         rhi = hi
>
>     res = where(reals < low,low,reals)
>     res = where(res > hi,hi,res)
>     return res.astype(a.typecode())
>
> I like the first one better.  It does *something* with an array no
> matter what its type is.  clip2() does to, but causes all code to run
> slower and is a lot more to write and test.  It is even more problematic
> on larger functions where the comparisons are buried within multiple
> computations.

I think you are drawing the wrong lesson here (but then I would). The second
problem does reveal a real deficiency in Numeric, but it's not the one you
think it is. It's the lack of support for generic access to the real and
imaginary parts of all arrays whether real or complex. This is easily coded,
but it's something that should really be in the Numeric core. I would write
the above functions as:

def clip_by_real(a, low, hi):
     a = asarray(a)
     real_a = real(a)
     res = where(real_a < low, low, a)
     res = where(real_a > hi , hi,  res)

This should have the same effect as clip on complex numbers where complex
comparisons are allowed. It has slightly different effects than clip2 which
always nukes the complex part of a number (probably not the correct
behaviour).

Of course real() and imag() don't currently exist, but they would be easy to
implement. In python a versions would be

def real(a):
     a = asarray(a)
     if a.typecode() in ['F','D']:
         return a.real
     else:
         return a

def imag(a):
     a = asarray(a)
     if a.typecode() in ['F','D']:
         return a.imag
     else:
         return np.zeros(a.shape, a.typecode())

Although it's possible efficiency could be improved, and if coded in C it
could certainly be improved.

-tim