\chapter{Ufuncs}
\label{cha:ufuncs}

\section{What are Ufuncs?}
\label{sec:what-are-ufuncs}

The operations on arrays that were mentioned in the previous section
(element-wise addition, multiplication, etc.) all share some features --- they
all follow similar rules for broadcasting, coercion and ``element-wise
operation''. Just as standard addition is available in Python through the add
function in the operator module, array operations are available through
callable objects as well. Thus, the following objects are available in the
numarray module:

\begin{table}[htbp]
   \centering
   \caption{Universal Functions, or ufuncs. The operators which invoke them
   when applied to arrays are indicated in parentheses. The entries in slanted
   typeface refer to unary ufuncs, while the others refer to binary ufuncs.} 
   \label{tab:ufuncs}
   \begin{tabular}{llll}
      add ($+$)         & subtract ($-$)   & multiply (*)       & divide ($/$) \\
      remainder (\%)    & power (**)       & \textsl{arccos}    & \textsl{arccosh} \\
      \textsl{arcsin}   & \textsl{arcsinh} & \textsl{arctan}    & \textsl{arctanh} \\
      \textsl{cos}      & \textsl{cosh}    & \textsl{tan}       & \textsl{tanh} \\
      \textsl{log10}    & \textsl{sin}     & \textsl{sinh}      & \textsl{sqrt} \\
      \textsl{absolute (abs)} & \textsl{fabs}    & \textsl{floor}     & \textsl{ceil} \\
      fmod              & \textsl{exp}     & \textsl{log}       & \textsl{conjugate} \\
      maximum           & minimum \\
      greater ($>$)     & greater\_equal ($>=$) & equal ($==$)  \\
      less ($<$)        & less\_equal ($<=$)  & not\_equal ($!=$) \\
      logical\_or       & logical\_xor     & logical\_not  & logical\_and \\
      bitwise\_or ($|$) & bitwise\_xor (\^{}) 
                        & bitwise\_not (\textasciitilde)  & bitwise\_and (\&)
      \\
      rshift ($>>$)       & lshift ($<<$)
   \end{tabular}
\end{table}

\remark{Add a remark here on how numarray does (or will) handle 'true'
and 'floor' division, which can be activated in Python 2.2 with
\samp{from __future__ import division}?.
Note: with 'true' division, \samp{1/2 == 0.5} and not \samp{0}.}

All of these ufuncs can be used as functions. For example, to use
\function{add}, which is a binary ufunc (i.e.\ it takes two arguments), one can
do either of:
\begin{verbatim}
>>> a = arange(10)
>>> print add(a,a)
[ 0  2  4  6  8 10 12 14 16 18]
>>> print a + a
[ 0  2  4  6  8 10 12 14 16 18]
\end{verbatim}
In other words, the \code{+} operator on arrays performs exactly the same thing
as the \function{add} ufunc when operated on arrays.  For a unary ufunc such as
\function{sin}, one can do, e.g.:
\begin{verbatim}
>>> a = arange(10)
>>> print sin(a)
[ 0.          0.84147096  0.90929741  0.14112    -0.7568025
      -0.95892429 -0.27941549  0.65698659  0.98935825  0.41211849]
\end{verbatim}
A unary ufunc returns an array with the same shape as its argument array, but
with each element replaced by the application of the function to that element
(\code{sin(0)=0}, \code{sin(1)=0.84147098}, etc.).

There are three additional features of ufuncs which make them different from
standard Python functions. They can operate on any Python sequence in addition
to arrays; they can take an ``output'' argument; they have methods which are
themselves callable with arrays and sequences. Each of these will be described
in turn.

Ufuncs can operate on any Python sequence. Ufuncs have so far been described as
callable objects which take either one or two arrays as arguments (depending on
whether they are unary or binary). In fact, any Python sequence which can be
the input to the \function{array} constructor can be used.  The return value
from ufuncs is always an array.  Thus:
\begin{verbatim}
>>> add([1,2,3,4], (1,2,3,4))
array([2, 4, 6, 8])
\end{verbatim}


\subsection{Ufuncs can take output arguments}
\label{sec:ufuncs-can-take}

In many computations with large sets of numbers, arrays are often used only
once. For example, a computation on a large set of numbers could involve the
following step
\begin{verbatim}
dataset = dataset * 1.20 
\end{verbatim}
This can also be written as the following using the Ufunc form:
\begin{verbatim}
dataset = multiply(dataset, 1.20)
\end{verbatim}
In both cases, a temporary array is created to store the results of the
computation before it is finally copied into \var{dataset}. It is
more efficient, both in terms of memory and computation time, to do an
``in-place'' operation.  This can be done by specifying an existing array as
the place to store the result of the ufunc. In this example, one can 
write:\footnote[1]{for Python-2.2.2 or later: `dataset *= 1.20' also works}
\begin{verbatim}
multiply(dataset, 1.20, dataset) 
\end{verbatim}
This is not a step to take lightly, however. For example, the ``big and slow''
version (\code{dataset = dataset * 1.20}) and the ``small and fast'' version
above will yield different results in at least one case:
\begin{itemize}
\item If the type of the target array is not that which would normally be
   computed, the operation will not coerce the array to the expected data type.
   (The result is done in the expected data type, but coerced back to the
   original array type.)
\item Example:
\begin{verbatim}
>>> a = arange(5, type=Float64)
>>> print a[::-1] * 1.2
[ 4.8  3.6  2.4  1.2  0. ]
>>> multiply(a[::-1], 1.2, a)
>>> a
array([ 4.8 ,  3.6 ,  2.4 ,  1.2,  0. ])
\end{verbatim}
\end{itemize}

The output array does not need to be the same variable as the input array. In
numarray, in contrast to Numeric, the output array may have any type (automatic
conversion is performed on the output).

\subsection{Ufuncs have special methods}
\label{sec:ufuncs-have-special-methods}


\begin{methoddesc}{reduce}{a, dim=0}
   If you don't know about the \function{reduce} command in Python, review
   section 5.1.3 of the Python Tutorial
   (\url{http://www.python.org/doc/current/tut/}). Briefly,
   \function{reduce} is most often used with two arguments, a callable object
   (such as a function), and a sequence. It calls the callable object with the
   first two elements of the sequence, then with the result of that operation
   and the third element, and so on, returning at the end the successive
   ``reduction'' of the specified callable object over the sequence elements.
   Similarly, the \method{reduce} method of ufuncs is called with a sequence as
   an argument, and performs the reduction of that ufunc on the sequence. As an
   example, adding all of the elements in a rank-1 array can be done with:
\begin{verbatim}
>>> a = array([1,2,3,4])
>>> print add.reduce(a)   # with Python's reduce, same as reduce(add, a)
10
\end{verbatim}
   When applied to arrays which are of rank greater than one, the reduction
   proceeds by default along the first axis:
\begin{verbatim}
>>> b = array([[1,2,3,4],[6,7,8,9]])
>>> print b
[[1 2 3 4]
 [6 7 8 9]]
>>> print add.reduce(b)
[ 7  9 11 13]
\end{verbatim}
   A different axis of reduction can be specified with a second integer
   argument:
\begin{verbatim}
>>> print b
[[1 2 3 4]
 [6 7 8 9]]
>>> print add.reduce(b, dim=1)
[10 30]
\end{verbatim}
\end{methoddesc}


\begin{methoddesc}{accumulate}{a}
   The \method{accumulate} ufunc method is simular to \method{reduce}, except
   that it returns an array containing the intermediate results of the
   reduction:
\begin{verbatim}
>>> a = arange(10)
>>> print a
[0 1 2 3 4 5 6 7 8 9]
>>> print add.accumulate(a)
[ 0  1  3  6 10 15 21 28 36 45] # 0, 0+1, 0+1+2, 0+1+2+3, ... 0+...+9
>>> print add.reduce(a) # same as add.accumulate(a)[-1] w/o side effects on a
45                                      
\end{verbatim}
\end{methoddesc}


\begin{methoddesc}{outer}{a, b}
   The third ufunc method is \method{outer}, which takes two arrays as
   arguments and returns the ``outer ufunc'' of the two arguments. Thus the
   \method{outer} method of the \function{multiply} ufunc, results in the outer
   product. The \method{outer} method is only supported for binary methods.
\begin{verbatim}
>>> print a
[0 1 2 3 4]
>>> print b
[0 1 2 3]
>>> print add.outer(a,b)
[[0 1 2 3]
 [1 2 3 4]
 [2 3 4 5]
 [3 4 5 6]
 [4 5 6 7]]
>>> print multiply.outer(b,a)
[[ 0  0  0  0  0]
 [ 0  1  2  3  4]
 [ 0  2  4  6  8]
 [ 0  3  6  9 12]]
>>> print power.outer(a,b)
[[ 1  0  0  0]
 [ 1  1  1  1]
 [ 1  2  4  8]
 [ 1  3  9 27]
 [ 1  4 16 64]]
\end{verbatim}
\end{methoddesc}


\begin{methoddesc}{reduceat}{}
   The reduceat method of Numeric has not been implemented in numarray.
\end{methoddesc}



\subsection{Ufuncs always return new arrays}
\label{sec:ufuncs-always-return}

Except when the output argument is used as described above, ufuncs always
return new arrays which do not share any data with the input arrays.


\section{Which are the Ufuncs?}
\label{sec:which-are-ufuncs}

Table \ref{tab:ufuncs} lists all the ufuncs. We will first discuss the
mathematical ufuncs, which perform operations very similar to the functions in
the \module{math} and \module{cmath} modules, albeit elementwise, on arrays.
These come in two forms, unary and binary:


\subsection{Unary Mathematical Ufuncs}
\label{sec:unary-math-ufuncs}

Unary ufuncs take only one argument.  The following ufuncs apply the
predictable functions on their single array arguments, one element at a time:
\function{arccos}, \function{arccosh}, \function{arcsin}, \function{arcsinh},
\function{arctan}, \function{arctanh}, \function{cos}, \function{cosh},
\function{exp}, \function{log}, \function{log10}, \function{sin},
\function{sinh}, \function{sqrt}, \function{tan}, \function{tanh},
\function{abs}, \function{fabs}, \function{floor}, \function{ceil},
\function{conjugate}.  As an example:
\begin{verbatim}
>>> print x
[0 1 2 3 4]
>>> print cos(x)
[ 1.          0.54030231 -0.41614684 -0.9899925  -0.65364362]
>>> print arccos(cos(x))
[ 0.          1.          2.          3.          2.28318531]
# not a bug, but wraparound: 2*pi%4 is 2.28318531
\end{verbatim}


\subsection{Binary Mathematical Ufuncs}
\label{sec:binary-math-ufuncs}

These ufuncs take two arrays as arguments, and perform the specified
mathematical operation on them, one pair of elements at a time: \function{add},
\function{subtract}, \function{multiply}, \function{divide},
\function{remainder}, \function{power}, \function{fmod}.


\subsection{Logical and bitwise ufuncs}
\label{sec:logical-ufuncs}

The ``logical'' ufuncs also perform their operations on arrays or numbers 
in elementwise fashion, just like the "mathematical" ones.  Two are special
(\function{maximum} and \function{miminum}) in that they return arrays with
entries taken from their input arrays:
\begin{verbatim}
>>> print x
[0 1 2 3 4]
>>> print y
[ 2.   2.5  3.   3.5  4. ]
>>> print maximum(x, y)
[ 2.   2.5  3.   3.5  4. ]
>>> print minimum(x, y)
[ 0.  1.  2.  3.  4.]
\end{verbatim}
The others logical ufuncs return arrays of 0's or 1's and of type Bool:
\function{logical_and}, \function{logical_or}, \function{logical_xor},
\function{logical_not}, 
These are fairly
self-explanatory, especially with the associated symbols from the standard
Python version of the same operations in Table \ref{tab:ufuncs} above. 
The bitwise ufuncs,
\function{bitwise_and}, \function{bitwise_or},
\function{bitwise_xor}, \function{bitwise_not},  
\function{lshift}, \function{rshift},  
on the other hand, only work with integer arrays (of any word size), and
will return integer arrays of the larger bit size of the two input arrays:
\begin{verbatim}
>>> x
array([7, 7, 0], type=Int8)
>>> y
array([4, 5, 6])
>>> x & y          # bitwise_and(x,y)
array([4, 5, 0])
>>> x | y          # bitwise_or(x,y)
array([7, 7, 6])   
>>> x ^ y          # bitwise_xor(x,y)
array([3, 2, 6]) 
>>> ~ x            # bitwise_not(x)
array([-8, -8, -1], type=Int8)

\end{verbatim}
We discussed finding contents of arrays based on arrays' indices by using slice.
Often, especially when dealing with the result of computations or data
analysis, one needs to ``pick out'' parts of matrices based on the content of
those matrices. For example, it might be useful to find out which elements of
an array are negative, and which are positive. The comparison ufuncs are
designed for such operation. Assume an array with various positive
and negative numbers in it (for the sake of the example we'll generate it from
scratch):
\begin{verbatim}
>>> print a
[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]
 [15 16 17 18 19]
 [20 21 22 23 24]]
>>> b = sin(a)
>>> print b
[[ 0.          0.84147098  0.90929743  0.14112001 -0.7568025 ]
 [-0.95892427 -0.2794155   0.6569866   0.98935825  0.41211849]
 [-0.54402111 -0.99999021 -0.53657292  0.42016704  0.99060736]
 [ 0.65028784 -0.28790332 -0.96139749 -0.75098725  0.14987721]
 [ 0.91294525  0.83665564 -0.00885131 -0.8462204  -0.90557836]]
>>> print greater(b, .3))
[[0 1 1 0 0]
 [0 0 1 1 1]
 [0 0 0 1 1]
 [1 0 0 0 0]
 [1 1 0 0 0]]
\end{verbatim}


\subsection{Comparisons}
\label{sec:comparisons}

The comparison functions \function{equal}, \function{not_equal},
\function{greater}, \function{greater_equal}, \function{less}, and
\function{less_equal} are invoked by the operators \code{==}, \code{!=},
\code{>}, \code{>=}, \code{<}, and \code{<=} respectively, but they can also be
called directly as functions. Continuing with the preceding example,
\begin{verbatim}
>>> print less_equal(b, 0)
[[1 0 0 0 1]
 [1 1 0 0 0]
 [1 1 1 0 0]
 [0 1 1 1 0]
 [0 0 1 1 1]]
\end{verbatim}
This last example has 1's where the corresponding elements are less than or
equal to 0, and 0's everywhere else.

The operators and the comparison functions are not exactly equivalent.  To
compare an array a with an object b, if b can be converted to an array, the
result of the comparison is returned. Otherwise, zero is returned. This means
that comparing a list and comparing an array can return quite different
answers. Since the functional forms such as equal will try to make arrays from
their arguments, using equal can result in a different result than using
\code{==}.
\begin{verbatim}
>>> a = array([1, 2, 3])
>>> b = [1, 2, 3]
>>> print a == 2
[0 1 0]
>>> print b == 2  
0          # (False since 2.3)
>>> print equal(a, 2)
[0 1 0]
>>> print equal(b, 2)
[0 1 0]
\end{verbatim}


\subsection{Ufunc shorthands}
\label{sec:ufunc-shorthands}

Numarray defines a few functions which correspond to popular ufunc methods:
for example, \function{add.reduce} is synonymous with the \function{sum}
utility function:
\begin{verbatim}
>>> a = arange(5)                       # [0 1 2 3 4]
>>> print sum(a)                        # 0 + 1 + 2 + 3 + 4
10
\end{verbatim}
Similarly, \function{cumsum} is equivalent to \function{add.accumulate} (for
``cumulative sum''), \function{product} to \function{multiply.reduce}, and
\function{cumproduct} to \function{multiply.accumulate}.  Additional
useful ``utility'' functions are \function{alltrue} and
\function{sometrue}, which are defined as \function{logical_and.reduce} and
\function{logical_or.reduce}, respectively:
\begin{verbatim}
>>> a = array([0,1,2,3,4])
>>> print greater(a,0)
[0 1 1 1 1]
>>> alltrue(greater(a,0))
0
>>> sometrue(greater(a,0))
1
\end{verbatim}

\section{Writing your own ufuncs!}
This section includes code from Examples/ufunc/setup.py in the numarray source
distribution to illustrate how to create a package which defines your own
universal functions.  Examples/ufunc is a working demonstration of what to do
and can also serve as a template to be modified.  Examples/ufunc in the
distribution source defines a particular directory structure which is assumed
by the following setup.py:

\begin{verbatim}
  ufunc
    |
    --- Lib
    |
    --- Src
\end{verbatim}

For a real extension, the top level directory should probably be named after
the extension rather than as it is named here, i.e. foo rather than ufunc.

\newpage
\begin{verbatim}
""" This setup demonstrates how to use the numarray code generation
facilities to define additional ufuncs at the user level.  Universal
functions apply operators or C-functions elementwise to all of the
elements of an array or pair of arrays. This setup generates code
defining three universal functions which are installed as the
numarray.foo package: Cos, bessel, and bar.

The ufuncs are used like this:

>>> import numarray as na
>>> import numarray.foo as foo
>>> a = na.array([1,2,3,4], type='Int32')
>>> foo.Cos(a)
array([ 0.54030228, -0.41614684, -0.9899925 , -0.65364361 ], type=Float32)

Note that since a is an Int32 array, Cos(a) first does blockwise
conversion of a to Float32 before calling the cosine library function.
"""

import distutils, os
from distutils.core import setup
from numarray.numarrayext import NumarrayExtension
from numarray.codegenerator import UfuncModule, make_stub

# ======================================================================
#
# Generate the C-code for the numarray.foo._foo extension:
# 

m = UfuncModule("_foo")            


# Inline the code for a couple C functions to be turned into ufuncs.
# This method lets you define your function here in the setup.py.  An
# alternate approach would be to link with a libarary or additional C
# module.
m.add_code("""

static double c_bessel(double x, double y)
{
      double x2=x*x, y2=y*y, diff=x2-y2;
      return diff*diff/(x2+y2);
}

static UInt8  c_bar(UInt32 x )
{
      return (x >> 24) & 0xFF;
}

""")

\end{verbatim}
\newpage
\begin{verbatim}

# Method add_ufunc() defines the name of a ufunc, it's corresponding C
# function which is applied element-wise to all elements of an array,
# The arity of the ufunc (unary or binary only), and the I/O type
# signatures which are directly supported by the ufunc.  Binary Ufunc
# "bessel" is implemented by the inline function "c_bessel" from the
# add_code() call above.

m.add_ufunc("bessel", "c_bessel", 2, [('Float32', 'Float32'),
                                      ('Float64', 'Float64')])

# Unary Ufunc "bar" only supports one builtin type signature:
# UInt32-->UInt8.  Other types are blockwise converted by the ufunc
# machinery to UInt32 before "c_bar" is called.

m.add_ufunc("bar", "c_bar", 1, [('UInt32','UInt8')])

# Unary Ufunc "cos" is implemented by the C standard library function
# "cos".  Given a Float32 array, it returns a Float32 array.  Given a
# Float64 array, it returns a Float64 array.  For other arrays,
# transparent conversion to one of the supported input types is performed
# by the ufunc machinery.

m.add_ufunc("cos", "cos", 1, [('Float32', 'Float32'),
                              ('Float64', 'Float64')])

# The generate() method writes out the complete extension module to the
# specified C file.
m.generate("Src/_foo.c")

# ======================================================================
# Create foo's __init__.py for defining UFuncs corresponding to CFuncs
# in _foo.  make_stub() emits boiler-plate which makes your extension
# a package.  The __init__ file imports all the public symbols from
# C extension _foo making them visible as numarray.foo.

make_stub("Lib/__init__", "_foo")  

# ======================================================================
#
# Standard distutils setup for the generated code.
#
setup(name = "foo",
      version = "0.1",
      maintainer = "Todd Miller",
      maintainer_email = "jmiller@stsci.edu",
      description = "foo universal functions for numarray",
      url = "http://www.stsci.edu/projects/foo/",
      packages = ["numarray.foo"],
      package_dir = { "numarray.foo":"Lib" },
      ext_modules = [ NumarrayExtension( 'numarray.foo._foo', ['Src/_foo.c'],
                                         # libraries = ['ttf'],
                                         # include_dirs = ['include'],
                                         # library_dirs = ['lib']
                                         )
                      ]
      )
\end{verbatim}


%% Local Variables:
%% mode: LaTeX
%% mode: auto-fill
%% fill-column: 79
%% indent-tabs-mode: nil
%% ispell-dictionary: "american"
%% reftex-fref-is-default: nil
%% TeX-auto-save: t
%% TeX-command-default: "pdfeLaTeX"
%% TeX-master: "numarray"
%% TeX-parse-self: t
%% End: