### 3.10 Declaring your own math functions on arrays ¶

There are four macros which make it easy to turn your own scalar functions into functions defined on arrays. They are:

    BZ_DECLARE_FUNCTION(f)                   // 1
    BZ_DECLARE_FUNCTION_RET(f,return_type)   // 2
    BZ_DECLARE_FUNCTION2(f)                  // 3
    BZ_DECLARE_FUNCTION2_RET(f,return_type)  // 4

Use version 1 when you have a function which takes one argument and returns a result of the same type. For example:

    #include <blitz/array.h>

    using namespace blitz;

    double myFunction(double x)
    {
        return 1.0 / (1 + x);
    }

    BZ_DECLARE_FUNCTION(myFunction)

    int main()
    {
        Array<double,2> A(4,4), B(4,4);  // ...
        B = myFunction(A);
    }

Use version 2 when you have a one argument function whose return type is different than the argument type, such as

    int g(double x);

Use version 3 for a function which takes two arguments and returns a result of the same type, such as:

    double g(double x, double y);

Use version 4 for a function of two arguments which returns a different type, such as:

    int g(double x, double y);

### 3.11 Tensor notation ¶

Blitz++ arrays support a tensor-like notation. Here’s an example of real-world tensor notation:

     ijk    ij k
    A    = B  C

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

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

    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;

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

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

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:

     ijk    ij k    jk i
    C    = A  x  - A  y

In Blitz++, this would be coded as:

    using namespace blitz::tensor;

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

This tensor expression can be visualized in the following way:

Examples of array indexing, subarrays, and slicing.

Here’s an example which computes an outer product of two one-dimensional arrays:

    #include <blitz/array.h>

    using namespace blitz;

    int main()
    {
        Array<float,1> x(4), y(4);
        Array<float,2> A(4,4);

        x = 1, 2, 3, 4;
        y = 1, 0, 0, 1;

        firstIndex i;
        secondIndex j;

        A = x(i) * y(j);

        cout << A << endl;

        return 0;
    }

And the output:

    (0,3) x (0,3)
    [ 1 0 0 1
      2 0 0 2
      3 0 0 3
      4 0 0 4 ]

Index placeholders can 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:

     ij    ik  kj
    C   = A   B

The repeated k index is interpreted as meaning

    c    = sum of (a   * b  ) over k
     ij             ik    kj

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

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

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

The `sum()` function is an example of an 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:

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

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

This is equivalent to the tensor notation

     ijkl    lj ki
    C     = A  B

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

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

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

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

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

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

The resulting implementation will look something like this:

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

The functions `cos` and `sin` will be invoked sixteen times each. It’s possible that a good optimizing compiler could hoist the `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 `cos(x)` and `sin(y)` first.

### 3.12 Array reductions ¶

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

- Reductions which transform an array into a scalar (for example, summing the elements). These are referred to as complete reductions.
- Reducing an N dimensional array (or array expression) to an N-1 dimensional array expression. These are called partial reductions.

### 3.13 Complete reductions ¶

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

    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

Here are the available complete reductions:

`sum()` ¶

  Summation (may be promoted to a higher-precision type)

`product()` ¶

  Product

`mean()` ¶

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

`min()` ¶

  Minimum value

`max()` ¶

  Maximum value

`minmax()` ¶

  Simultaneous minimum and maximum value (returns a value of type MinMaxValue<T>)

`minIndex()` ¶

  Index of the minimum value (`TinyVector<int,N_rank>`)

`maxIndex()` ¶

  Index of the maximum value (`TinyVector<int,N_rank>`)

`count()` ¶

  Counts the number of times the expression is logical true (`int`)

`any()` ¶

  True if the expression is true anywhere (`bool`)

`all()` ¶

  True if the expression is true everywhere (`bool`)

Caution: `minIndex()` and `maxIndex()` return TinyVectors, even when the rank of the array (or array expression) is 1.

Reductions can be combined with `where` expressions (where statements) to reduce over some part of an array. For example, `sum(where(A > 0, A, 0))` sums only the positive elements in an array.

### 3.14 Partial Reductions ¶

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

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

    rs = sum(A, j);

The reduction `sum()` takes two arguments:

- The first argument is an array or array expression.
- The second argument is an index placeholder indicating the dimension over which the reduction is to occur.

Reductions have an important restriction: It is currently only possible to reduce over the 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:

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

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

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:

    B = sum(A(i,k,j), k);

Writing `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:

`sum()`

  Summation

`product()`

  Product

`mean()`

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

`min()`

  Minimum value

`max()`

  Maximum value

`minIndex()`

  Index of the minimum value (int)

`maxIndex()`

  Index of the maximum value (int)

`count()`

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

`any()`

  True if the expression is true anywhere (bool)

`all()`

  True if the expression is true everywhere (bool)

`first()`

  First index at which the expression is logical true (int); if the expression is logical true nowhere, then `tiny(int())` (INT_MIN) is returned.

`last()`

  Last index at which the expression is logical true (int); if the expression is logical true nowhere, then `huge(int())` (INT_MAX) is returned.

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

To illustrate, here’s an 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(4);
    firstIndex i;
    secondIndex j;

    z = sum(A(j,i), j);

The array `z` now contains the sum of `A` along each column:

    [ 10    5     12    7 ]

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

    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 ]

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

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

    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);

Note that this is not allowed:

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

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

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