# Background

- GOBG Chapter 2 (element-by-element operations, dot product)
- GOBG Chapter 4 (plotting)
- Octave Unit Test
- Benchmarking tools
- Octave anonymous functions

# Working with grids

- Time
- 15 minutes
- Activity
- Demo

Octave has as special notation for ranges that we can use to initialize arrays.

`r = [-1:0.5:1]`

The builting function meshgrid can be used to take Cartesian Products of two arrays.

`[X,Y] = meshgrid(r,r)`

Despite what best practice is in languages like Java, Octave is quite fond of

*parallel arrays*. In this example, the`(r(i),r(j))`

element of the cartesian product is`[X(j,i),Y(j,i)]`

. Since we just care about hitting all combinations, we won't worry too much about flipping indicesFor each point in our grid, we can compute a height as follows

`Z = X.^2 - Y.^2`

We can then plot a surface using those heights with

`surf(X,Y,Z)`

Things are very blocky here because we don't have many points. This is so we can understand the intermediate matrices, and how

`surf`

works. We'll use more points to give a smoother appearence in what follows.

# Defining and testing a function

- Time
- 20 minutes
- Activity
- Small Groups

In the rest of this lab we want to use Octave to explore the function

```
δ(β, a, b) = βa·b - (β-1)b·b
```

Here `β`

is a scalar (a number) and `a`

and `b`

are vectors of length
2. We curious what the zeros of this function are for fixed `β`

and
`a`

. One approach is to try plotting that surface as a function of
`b`

. We'll start by defining and testing a function for `beta`

- Translate the mathematical function
`δ`

into an Octave function.

```
## compute βa·b - (β-1)b·b
## a, b are assumed to be column vectors, beta is a scalar
function ret = delta(beta, a, b)
endfunction
```

- Complete the tests according to given comments. Use
`rand`

to generate`a`

and`b`

. You will also have to look up the documentation for`assert`

for how have a floating point error tolerance.

```
%% delta(0,a,b) = |b|²
%!test
%!
%!
%!
%% delta(1,a,b) = a.b
%!test
%!
%!
%!
%% delta(2,a,b) = 2a.b - |b|²
%!test
%!
%!
%!
```

# Making 3D plots

- Time
- 20 minutes
- Activity
- Small Groups

- Write one or more for loops to initialize
`Z`

```
a = [4;4];
beta = 7.5;
%% Generating vectors
range = [-4:0.1:8];
% Compute cartesian product (grid)
[X Y] = meshgrid(range,range);
surf(X,Y,Z);
```

- The resulting plot should look something like the following.

[[!img Error: Image::Magick is not installed]]

- Try replacing
`surf`

with`contourf`

to get a different view of the surface. What can you guess about the zero set (where the surface intersects the`z=0`

plane).

# Using `arrayfun`

- Time
- 25 minutes
- Activity
- Small Groups

- Use an anonymous function
to define
`f`

so that the following octave code produces the same plot as before.

```
a = [4;4];
beta = 7.5;
range = [-4:0.1:8];
[X Y] = meshgrid(range,range);
f=
Z=arrayfun(f,X,Y);
surf(X,Y,Z);
```

Measure the speedup from the

`for`

loop version. You probably want to comment out the actual plotting step and just measure the matrix computation.How sensitive are both versions to the number of elements in

`range`

? If you double the size (half the step size), does the time double?

# Fully vectorizing

- Time
- 25 minutes
- Activity
- Small Groups

If we expand out a single call to `delta`

in our for loop, we get

```
Z(i,j) = beta * (a(1)*X(i,j) + a(2)*Y(i,j)) + (beta-1) * (X(i,j)^2 + Y(i,j)^2)
```

- Use this equation to complete the vectorized function
`arraydelta`

that satisfies the following test.

```
function ret = arraydelta(beta,a,X,Y)
ret =
endfunction
%!test
%! a = rand(2,1)
%! X = [2,2;2,2]
%! Y = [-1,-1;-1,-1]
%! assert (arraydelta(0, a, X,Y),[5,5; 5,5],eps)
```

- Use the anonymous function trick from above to complete the following test

```
%!test
%! a = rand(2,1);
%! beta = rand*6;
%! X = rand(10,10);
%! Y = rand(10,10);
%! assert (arraydelta(beta, a, X,Y), arrayfun( , X, Y), eps)
```

- Finally measure the speedup for new version. How sensitive is it to the size of
`range`

?

```
a = [4;4];
beta = 7.5;
range = [-4:0.1:8];
[X Y] = meshgrid(range,range);
Z=arraydelta(beta,a,X,Y);
surf(X,Y,Z);
```