arrays-expr.texi

A C++ class library for scientific computing

@end ifnottex

@math{A} is a rank 3 tensor (a three dimensional array), @math{B} is a rank
2 tensor (a two dimensional array), and @math{C} is a rank 1 tensor (a one
dimensional array).  The above expression sets 
@code{A(i,j,k) = B(i,j) * C(k)}.

To implement this product using Blitz++, we'll need the arrays and some
index placeholders:

@cindex index placeholders used for tensor notation

@example
Array<float,3> A(4,4,4);
Array<float,2> B(4,4);
Array<float,1> C(4);

firstIndex i;    // Alternately, could just say
secondIndex j;   // using namespace blitz::tensor;
thirdIndex k;
@end example

Here's the Blitz++ code which is equivalent to the tensor expression:

@example
A = B(i,j) * C(k);
@end example

The index placeholder arguments tell an array how to map its dimensions onto
the dimensions of the destination array.  For example, here's some
real-world tensor notation:
@tex
$$ C^{ijk} = A^{ij} x^{k} - A^{jk} y^{i} $$
@end tex
@html
<pre>
 ijk    ij k    jk i
C    = A  x  - A  y
</pre>
@end html
@ifnottex
@ifnothtml
@example
 ijk    ij k    jk i
C    = A  x  - A  y
@end example
@end ifnothtml
@end ifnottex

In Blitz++, this would be coded as:

@example
using namespace blitz::tensor;

C = A(i,j) * x(k) - A(j,k) * y(i);
@end example

This tensor expression can be visualized in the following way:

@center @image{tensor1}
@center Examples of array indexing, subarrays, and slicing.

Here's an example which computes an outer product of two one-dimensional
arrays:
@cindex outer product
@cindex kronecker product
@cindex tensor product

@smallexample
@include examples/outer.texi
@end smallexample

And the output:

@smallexample
@include examples/outer.out
@end smallexample

Index placeholders can @emph{not} be used on the left-hand side of an
expression.  If you need to reorder the indices, you must do this on the
right-hand side.

In real-world tensor notation, repeated indices imply a contraction (or
summation).  For example, this tensor expression computes a matrix-matrix
product:
@tex
$$ C^{ij} = A^{ik} B^{kj} $$
@end tex
@html
<pre>
 ij    ik  kj
C   = A   B
</pre>
@end html
@ifnottex
@ifnothtml
@example
 ij    ik  kj
C   = A   B
@end example
@end ifnothtml
@end ifnottex

The repeated k index is interpreted as meaning
@tex
$$ c_{ij} = \sum_{k} a_{ik} b_{kj} $$
@end tex
@html
<pre>
c    = sum of {a   * b  } over k
 ij             ik    kj
</pre>
@end html
@ifnottex
@ifnothtml
@example
c    = sum of @{a   * b  @} over k
 ij              ik    kj
@end example
@end ifnothtml
@end ifnottex

@cindex contraction
@cindex tensor contraction

In Blitz++, repeated indices do @emph{not} imply contraction.  If you want
to contract (sum along) an index, you must use the @code{sum()} function:

@example
Array<float,2> A, B, C;   // ...
firstIndex i;
secondIndex j;
thirdIndex k;

C = sum(A(i,k) * B(k,j), k);
@end example

The @code{sum()} function is an example of an @emph{array reduction},
described in the next section.

Index placeholders can be used in any order in an expression.  This example
computes a kronecker product of a pair of two-dimensional arrays, and
permutes the indices along the way:

@example
Array<float,2> A, B;   // ...
Array<float,4> C;      // ...
fourthIndex l;

C = A(l,j) * B(k,i);
@end example

This is equivalent to the tensor notation
@tex
$$ C^{ijkl} = A^{lj} B^{ki} $$
@end tex
@html
<pre>
 ijkl    lj ki
C     = A  B
 </pre>
@end html
@ifnottex
@ifnothtml
@example
 ijkl    lj ki
C     = A  B
@end example
@end ifnothtml
@end ifnottex

Tensor-like notation can be mixed with other array notations:

@example
Array<float,2> A, B;  // ...
Array<double,4> C;    // ...

C = cos(A(l,j)) * sin(B(k,i)) + 1./(i+j+k+l);
@end example

@cindex tensor notation efficiency issues
An important efficiency note about tensor-like notation: the right-hand side
of an expression is @emph{completely evaluated} for @emph{every} element in
the destination array.  For example, in this code:

@example
Array<float,1> x(4), y(4);
Array<float,2> A(4,4):

A = cos(x(i)) * sin(y(j));
@end example

The resulting implementation will look something like this:

@example
for (int n=0; n < 4; ++n)
  for (int m=0; m < 4; ++m)
    A(n,m) = cos(x(n)) * sin(y(m));
@end example

The functions @code{cos} and @code{sin} will be invoked sixteen times each.
It's possible that a good optimizing compiler could hoist the @code{cos}
evaluation out of the inner loop, but don't hold your breath -- there's a
lot of complicated machinery behind the scenes to handle tensor notation,
and most optimizing compilers are easily confused.  In a situation like the
above, you are probably best off manually creating temporaries for
@code{cos(x)} and @code{sin(y)} first.

@section Array reductions
@cindex Array reductions
@cindex reductions

Currently, Blitz++ arrays support two forms of reduction:

@itemize @bullet

@item      Reductions which transform an array into a scalar (for example,
summing the elements).  These are referred to as @strong{complete
reductions}.

@item      Reducing an N dimensional array (or array expression) to an N-1
dimensional array expression.  These are called @strong{partial reductions}.

@end itemize

@cindex Array reductions complete
@cindex complete reductions
@cindex reductions complete

@section Complete reductions

Complete reductions transform an array (or array expression) into 
a scalar.  Here are some examples:

@example 
Array<float,2> A(3,3);
A = 0, 1, 2,
    3, 4, 5,
    6, 7, 8;
cout << sum(A) << endl          // 36
     << min(A) << endl          // 0
     << count(A >= 4) << endl;  // 5
@end example

Here are the available complete reductions:

@table @code
@item sum()
@cindex @code{sum()} reduction
Summation (may be promoted to a higher-precision type)

@item product()
@cindex @code{product()} reduction
Product 

@item mean()
@cindex @code{mean()} reduction
Arithmetic mean (promoted to floating-point type if necessary) 

@item min()
@cindex @code{min()} reduction
Minimum value 

@item max()
@cindex @code{max()} reduction
Maximum value 

@item minIndex()
@cindex @code{minIndex()} reduction
Index of the minimum value (@code{TinyVector<int,N_rank>})

@item maxIndex()
@cindex @code{maxIndex()} reduction
Index of the maximum value (@code{TinyVector<int,N_rank>})

@item count()
@cindex @code{count()} reduction
Counts the number of times the expression is logical true (@code{int})

@item any()
@cindex @code{any()} reduction
True if the expression is true anywhere (@code{bool})

@item all()
@cindex @code{all()} reduction
True if the expression is true everywhere (@code{bool})
@end table

@strong{Note:} @code{minIndex()} and @code{maxIndex()} return TinyVectors, 
even when the rank of the array (or array expression) is 1.

Reductions can be combined with @code{where} expressions (@ref{Where expr})
to reduce over some part of an array.  For example, @code{sum(where(A > 0,
A, 0))} sums only the positive elements in an array.

@section Partial Reductions

@cindex Array reductions partial
@cindex partial reductions
@cindex reductions partial

Here's an example which computes the sum of each row of a two-dimensional
array:

@example
Array<float,2> A;    // ...
Array<float,1> rs;   // ...
firstIndex i;
secondIndex j;

rs = sum(A, j);
@end example

The reduction @code{sum()} takes two arguments:

@itemize @bullet

@item     The first argument is an array or array expression.

@item     The second argument is an index placeholder indicating the
dimension over which the reduction is to occur.  

@end itemize

Reductions have an @strong{important restriction}: It is currently only
possible to reduce over the @emph{last} dimension of an array or array
expression.  Reducing a dimension other than the last would require Blitz++
to reorder the dimensions to fill the hole left behind.  For example, in
order for this reduction to work:

@example
Array<float,3> A;   // ...
Array<float,2> B;   // ...
secondIndex j;

// Reduce over dimension 2 of a 3-D array?
B = sum(A, j);
@end example

Blitz++ would have to remap the dimensions so that the third dimension
became the second.  It's not currently smart enough to do this.

However, there is a simple workaround which solves some of the problems
created by this limitation: you can do the reordering manually, prior to the
reduction:

@example
B = sum(A(i,k,j), k);
@end example

Writing @code{A(i,k,j)} interchanges the second and third dimensions,
permitting you to reduce over the second dimension.  Here's a list of the
reduction operations currently supported:

@table @code
@item sum()    
Summation

@item product() 
Product 

@item mean()   
Arithmetic mean (promoted to floating-point type if necessary)

@item min()    
Minimum value

@item max()    
Maximum value

@item minIndex() 
Index of the minimum value (int)

@item maxIndex() 
Index of the maximum value (int)

@item count()    
Counts the number of times the expression is logical true (int)

@item any()      
True if the expression is true anywhere (bool)

@item all()      
True if the expression is true everywhere (bool)

@item first()   
First index at which the expression is logical true (int); if the expression
is logical true nowhere, then @code{tiny(int())} (INT_MIN) is returned.

@item last()     
Last index at which the expression is logical true (int); if the expression
is logical true nowhere, then @code{huge(int())} (INT_MAX) is returned.  
@end table

The reductions @code{any()}, @code{all()}, and @code{first()} have
short-circuit semantics: the reduction will halt as soon as the answer is
known.  For example, if you use @code{any()}, scanning of the expression
will stop as soon as the first true value is encountered.

To illustrate, here's an example:

@example
Array<int, 2> A(4,4);

A =  3,   8,   0,   1,
     1,  -1,   9,   3,
     2,  -5,  -1,   1,
     4,   3,   4,   2;

Array<float, 1> z;
firstIndex i;
secondIndex j;

z = sum(A(j,i), j);
@end example

The array @code{z} now contains the sum of @code{A} along each column:

@example
[ 10    5     12    7 ]
@end example

This table shows what the result stored in @code{z} would be if
@code{sum()} were replaced with other reductions:

@example
sum                     [         10         5        12         7 ]
mean                    [        2.5      1.25         3      1.75 ]
min                     [          1        -5        -1         1 ]
minIndex                [          1         2         2         0 ]
max                     [          4         8         9         3 ]
maxIndex                [          3         0         1         1 ]
first((A < 0), j)       [ -2147483648        1         2 -2147483648 ]
product                 [         24       120         0         6 ]
count((A(j,i) > 0), j)  [          4         2         2         4 ]
any(abs(A(j,i)) > 4, j) [          0         1         1         0 ]
all(A(j,i) > 0, j)      [          1         0         0         1 ]
@end example

Note: the odd numbers for first() are @code{tiny(int())} i.e.@: the smallest
number representable by an int.  The exact value is machine-dependent.

@cindex Array reductions chaining
@cindex partial reductions chaining
@cindex reductions chaining

The result of a reduction is an array expression, so reductions
can be used as operands in an array expression:

@example
Array<int,3> A;
Array<int,2> B;
Array<int,1> C;   // ...

secondIndex j;
thirdIndex k;

B = sqrt(sum(sqr(A), k));

// Do two reductions in a row
C = sum(sum(A, k), j);
@end example

Note that this is not allowed:

@example
Array<int,2> A;
firstIndex i;
secondIndex j;

// Completely sum the array?
int result = sum(sum(A, j), i);
@end example

You cannot reduce an array to zero dimensions!  Instead, use one of the
global functions described in the previous section.


@node Where expr
@section where statements
@cindex @code{where} statements
@cindex functional if (@code{where})
@cindex @code{if} (@code{where})

Blitz++ provides the @code{where} function as an array expression version of the
@code{( ? : )} operator.  The syntax is:

@example
where(array-expr1, array-expr2, array-expr3)
@end example

Wherever @code{array-expr1} is true, @code{array-expr2} is returned.  Where
@code{array-expr1} is false, @code{array-expr3} is returned.  For example,
suppose we wanted to sum the squares of only the positive elements of an
array.  This can be implemented using a where function:

@example
double posSquareSum = sum(where(A > 0, pow2(A), 0));
@end example