# Category Archives: Scientific Computing

## N-dim Array Broadcasting Efficiency in ILNumerics 4.10

Next to other great improvements in version 4.10 of ILNumerics Ultimate VS, it is especially one new feature which requires some attention: general broadcasting for n-dimensional ILArrays.

Broadcasting as a concept today is found in many popular mathematical prototyping systems. The most direct correspondence probably exists in the numpy package. Matlab and Octave offer similar functionality by means of the bsxfun function.

The term ‘broadcasting’ refers to a binary operator which is able to apply an elementwise operation to the elements of two n-dimensional arrays. The ‘operation’ is often as simple as a straight addition (ILMath.add, ILMath.divide, a.s.o.). However, what is special about broadcasting is that it allows the operation even for the case where both arrays involved do not have the same number of elements.

## Broadcasting in ILNumerics prior Version 4.10

In ILNumerics, broadcasting is available for long already. But prior version 4.10 it was limited to scalars operating on n-dim arrays and vectors operating on matrices. Therefore, we had used the term ‘vector expansion’ instead of broadcasting. Obviously, broadcasting can be seen as a generalization of vector expansion.

Let’s visualize the concept by considering the following matrix A:

A

1  5   9  13  17
2  6  10  14  18
3  7  11  15  19
4  8  12  16  20


Matrix A might represent 5 data points of a 4 dimensional dataset as columns. One common requirement is to apply a certain operation to all datapoints in a similar way. In order to, let’s say, scale/weight the components of each dimension by a certain factor, one would multiply each datapoint with a vector of length 4.


ILArray<double> V = new[] { 0.5, 3.0, 0.5, 1.0 };

0.5
3.0
0.5
1.0


The traditional way of implementing this operation would be to expand the scaling vector by replicating it from a single column to a matrix matching the size of A.

VExp = repmat(V, 1, 5);

0.5  0.5  0.5  0.5  0.5
3.0  3.0  3.0  3.0  3.0
0.5  0.5  0.5  0.5  0.5
1.0  1.0  1.0  1.0  1.0


Afterwards, the result can be operated with A elementwise in the common way.

ILArray<double> Result = VExp * A;

0.5   2.5   4.5   6.5   8.5
6.0  18.0  30.0  42.0  54.0
1.5   3.5   5.5   7.5   9.5
4.0   8.0  12.0  16.0  20.0


The problem with the above approach is that the vector data need to be expanded first. There is little advantage in doing so: a lot of new memory is being used up in order to store completely redundant data. We all know that memory is the biggest bottleneck today. We should prevent from lots of memory allocations whenever possible. This is where vector expansion comes into play. In ILNumerics, for long, one can prevent from the additional replication step and operate the vector on the matrix directly. Internally, the operation is implemented in a very efficient way, without replicating any data, without allocating new memory.

ILArray<double> Result = V * A;

0.5   2.5   4.5   6.5   8.5
6.0  18.0  30.0  42.0  54.0
1.5   3.5   5.5   7.5   9.5
4.0   8.0  12.0  16.0  20.0

## Generalizing for n-Dimensions

Representing data as matrices is very popular in scientific computing. However, if the data are stored into arrays of other shapes, having more than two dimensions, one had to fall back to repmatting in order for the binary operation to succeed. This nuissance has been removed in version 4.10.

Now it is possible to apply broadcasting to two arrays of any matching shape – without the need for using repmat. In order for two arrays to ‘match‘ in the binary operation, the following rules must be fullfilled:

1. All corresponding dimensions of both arrays must match.
2. In order for two  corresponding dimensions to match,
• both dimensions must be of the same length, or
• one of the dimensions must be of length 1.

An example of two matching arrays would be a vector running along the 3rd dimension and a 3 dimensional array:

In the above image the vector (green) has the same length as the corresponding dimension of the 3D array (gray). The size of the vector is [1 x 1 x 6]. The size of the 3D array is [4 x 5 x 6]. Hence, any dimension of both, the vector and the 3D array ‘match’ in terms of broadcasting. A broadcasting operation for both, the vector and the array would give the same result as if the vector would be replicated along the 1st and the 2nd dimensions. The first element will serve all elements in the first 4 x 5 slice in the 1-2 plane. This slice is marked red in the next image: Note that all red elements here derive from the same value – from the first element of the green vector.  The same is true for all other vector elements: they fill corresponding slices on the 3D array along the 3rd dimension.

Slowly, a huge performance advantage of broadcasting becomes clear: the amount of memory saved explodes when more, longer dimensions are involved.

## Special Case: Broadcasting on Vectors

In the most general case and if broadcasting is blindly applied, the following special case potentially causes issues. Consider two vectors, one row vector and one column vector being provided as input parameters to a binary operation. In ILNumerics, every array carries at least two dimensions. A column vector of length 4 is considered an array of size [4 x 1]. A row vector of length 5 is considered an array of size [1 x 5]. In fact, any two vectors match according to the general  broadcasting rules.

As a consequence operating a row vector [1 x 5] with a column vector [4 x 1] results in a matrix [4 x 5]. The row vector is getting ‘replicated’ (again, without really executing the replication) four times along the 1st dimension, and the column vector 5 times along the rows.

array(new[] {1.0,2.0,3.0,4.0,5.0}, 1, 5) +array(new[] {1.0,2.0,3.0,4.0}, 4, 1)

<Double> [4,5]
[0]:          2          3          4          5          6
[1]:          3          4          5          6          7
[2]:          4          5          6          7          8
[3]:          5          6          7          8          9


Note, in order for the above code example to work, one needs to apply a certain switch:

Settings.BroadcastCompatibilityMode = false;


The reason is that in the traditional version of ILNumerics (just like in Matlab and Octave) the above code would simply not execute but throw an exception instead. Originally, binary operations on vectors would ignore the fact that vectors are matrices and only take the length of the vectors into account, operating on the corresponding elements if the length of both vectors do match. Now, in order to keep compatibility for existing applications, we kept the former behavior.

The new switch ‘Settings.BroadcastCompatibilityMode’ by default is set to ‘true’. This will cause the Computing Engine to throw an exception when two vectors of inequal length are provided to binary operators. Applying vectors of the same length (regardless of their orientation) will result in a vector of the same length.

If the ‘Settings.BroadcastCompatibilityMode’ switch is set to ‘false’ then general broadcasting is applied in all cases according to the above rules – even on vectors. For the earlier vector example this leads to the resulting matrix as shown above: operating a row on a column vector expands both vectors and gives a matrix of corresponding size.

Further reading: binary operators, online documentation

## How did you do that? – The Computation

Wrapping up our series of ILNumerics and science we take a closer look at computation and science. Easily code your science into a robust, reliable program with high-speed performance.

## ILNumerics for Scientists – Going 3D

### Recap

Last time I started with one of the easiest problems in quantum mechanics: the particle in a box. This time I’ll add 1 dimension and we’ll see a particle in a 2D box. To visualize its wave function and density we need 3D surface plots.

### 2D Box

This time we have a particle that is confined in a 2D box. The potential within the box is zero and outside the box infinity. Again the solution is well-known and can be found on Wikipedia. This time the state of the wave function is determined by two numbers. These are typically called quantum numbers and refer to the X and the Y direction, respectively.

The absolute size of the box doesn’t really matter and we didn’t worry about it in the 1D case. However, the relative size of the length and the width make a difference. The solution to our problem reads

$\Psi_{n,k}(x,y) = \sqrt{\frac{4}{L_x L_y}} \cdot \sin(n \cdot \pi \cdot x / L_x) \cdot \sin(k \cdot \pi \cdot y / L_y)$

### The Math

Very similar to the 1D case I quickly coded the wave function and the density for further plotting. I had to make sure that the arrays are fit for 3D plotting, so the code looks a little bit different compared to last post’s

     public static ILArray<double> CalcWF(int EVXID, int EVYID, double LX, double LY, int MeshSize)
{
ILArray<double> X = linspace<double>(0, LX, MeshSize);
ILArray<double> Y = linspace<double>(0, LY, MeshSize);

ILArray<double> Y2d = 1;
ILArray<double> X2d = meshgrid(X, Y, Y2d);

ILArray<double> Z = sqrt(4.0 / LX / LY) * sin(EVXID * pi * X2d / LX) * sin(EVYID * pi * Y2d / LY);

return Z.Concat(X2d,2).Concat(Y2d,2);
}


Again, this took me like 10 minutes and I was done.

### The Visualization

This time the user can choose the quantum numbers for X and Y direction, the ratio between the length and the width of the box and also the number of mesh points along each axis for plotting. This makes the visualization panel a little bit more involved. Nevertheless, it’s still rather simple and easy to use. This time it took me only 45 minutes – I guess I learned a lot from last time.

### The result

Here is the result of my little program. You can click and play with it. If you’re interested, you can download the Particle2DBox source code. Have fun!

This is a screenshot of the application. I chose the second quantum number along the x axis and the fourth quantum number along the y axis. The box is twice as long in y direction as it is in x direction. The mesh size is 100 in each direction. On the left hand side you see the wave function and on the right hand side the probability density.

## Directions to the ILNumerics Optimization Toolbox

As of yesterday the ILNumerics Optimization Toolbox is out and online! It’s been quite a challenge to bring everything together: some of the best algorithms, the convenience you as a user of ILNumerics expect and deserve, and the high performance requirements ILNumerics sets the scale on for. We believe that all these goals could be achieved quite greatly.

## ILNumerics for Scientists – An easy start

### Motivation

I’ve been working as a scientist at universities for 10 years before deciding to go into industry. The one thing I hated most was coding. At the end of the day coding for scientists is like running for a football player. Obviously, you need it but it’s not what you’re here for.

I really dreaded the coding and the debugging. So much precious time for something that was so clear on paper and I just wanted the solution of my equations to see whether my idea made sense or not. More often than not scientists find that their idea was not so great and now they had spent so much time coding just to find out that the idea didn’t work. Continue reading ILNumerics for Scientists – An easy start

## C# for 3D visualizations and Plotting in .NET

2D and 3D Visualizations are an important feature for a wide range of domains: both software developers and scientists often need convenient visualization facilities to create interactive scenes and to make data visible. The ILNumerics math library brings powerful visualization features to C# and .NET: ILView, the ILNumerics Scene Graph API and its plotting engine. We’d like to give an overview over our latest achievements.

## ILView: a simple way to create interactive 3d visualizations

We have created ILView as an extension to our interactive web component: It allows you to simply try out ILNumerics’ 2d and 3d visualization features by chosing the output format .exe in our visualization examples. But that’s not all: ILView is also a general REPL for the evaluation of computational expressions using C# language. ILView is Open Source – find it on GitHub!

## ILNumerics Scene Graph: realize complex visualizations in .NET

The ILNumeric’s scene graph is the core of ILNumerics’ visualization engine. No matter if you want to create complex interactive 3D visualizations, or if you aim at enhancing and re-configuring existing scenes in .NET: The ILNumerics scene graph offers a convenient way to realize stunning graphics with C#. It uses OpenGL, GDI, and it’s possible to export scenes into vector and pixel graphics.

## Scientific Plotting: visualize your data using C#

With ILNumerics’ visualization capabilities, C# becomes the language of choice for scientists, engineers and developers who need to visualize data: Our plotting API and different kinds of plotting types (contour plots, surface plots etc.) make easy work of creating beautiful scientific visualizations.