---
engine: julia
---

```{julia}
#| error: false
#| echo: false
#| output: false
using InteractiveUtils
import QuartoNotebookWorker
Base.stdout = QuartoNotebookWorker.with_context(stdout)
myactive_module() = Main.Notebook
Base.active_module() = myactive_module()
```

# Vectors, Matrices, Arrays

## General 

Let us now turn to the probably most important containers for numerical mathematics: 

- Vectors `Vector{T}` 
- Matrices `Matrix{T}` with two indices  
- N-dimensional arrays with N indices `Array{T,N}`

In fact, `Vector{T}` is an alias for `Array{T,1}` and `Matrix{T}` is an alias for `Array{T,2}`.


```{julia}
Vector{Float64} === Array{Float64,1} && Matrix{Float64} ===  Array{Float64,2}
```

When created by an explicit element list, the 'greatest common type' for the type parameter `T` is determined.

```{julia}
v = [33, "33", 1.2]
```
If `T` is a numeric type, the elements are converted to this type.

```{julia}
v = [3//7, 4, 2im]
``` 

The following functions provide information about an array:

- `length(A)` --- number of elements
- `eltype(A)`  --- type of elements
- `ndims(A)`   --- number of dimensions (indices)
- `size(A)`   ---   tuple with the dimensions of the array  

```{julia}
v1 = [12, 13, 15]

m1 = [ 1    2.5
       6    -3 ]

for f ∈ (length, eltype, ndims, size)
    println("$(f)(v) = $(f(v)),          $(f)(m1) = $(f(m1))")
end
```

- The strength of the 'classical' array for scientific computing lies in the fact that it is simply a contiguous memory segment in which components of the same length (e.g. 64 bits) are stored sequentially. This makes the memory requirement minimal and the access speed to a component, both when reading and when modifying, maximal. The location of the component `v[i]` can be calculated immediately from `i`. 
- Julia's `Array{T,N}` (and therefore vectors and matrices) is implemented in this way for the usual numeric types `T`. The elements are stored *unboxed*. In contrast, for example, a `Vector{Any}` is implemented as a list of addresses of objects *(boxed)* and not as a list of the objects themselves. 
- Julia's `Array{T,N}` stores its elements directly *(unboxed)*, when `isbitstype(T) == true`.

```{julia}
isbitstype(Float64), 
isbitstype(Complex{Rational{Int64}}), 
isbitstype(String)
```

## Vectors

### List-like Functions 

- `push!(vector, items...)` --- appends elements at the end of the vector
- `pushfirst!(vector, items...)` --- prepends elements at the beginning of the vector
- `pop!(vector)` --- removes the last element and returns it as a result,
- `popfirst!(vector)` --- removes the first element and returns it  

```{julia}
v = Float64[]           # empty Vector{Float64}

push!(v, 3, 7)
pushfirst!(v, 1)

a = pop!(v)
println("a= $a")

push!(v, 17)
```

A `push!()` can be very expensive, as new memory may need to be allocated and then the entire existing vector needs to be copied. Julia optimizes the memory management. In such a case, memory is allocated in advance, so that further `push!`s are very fast and one 'almost achieves O(1) speed'.   

However, one should avoid operations like `push!()` or `resize()` in time-critical code and with very large arrays.

### Further Constructors

One can create vectors with a given length and type uninitialized. This is fastest, the elements are random bit patterns.
```{julia}
# fixed length 1000, uninitialized

v = Vector{Float64}(undef, 1000)
v[345]
```


- `zeros(n)` creates a `Vector{Float64}` of length `n` and initializes with zero.  

```{julia}
v = zeros(7) 
```


- `zeros(T,n)` creates a zero vector of type `T`.

```{julia}
v=zeros(Int, 4)
```

- `fill(x, n)` creates a `Vector{typeof(x)}` of length `n` and fills with `x`.

```{julia}
v = fill(sqrt(2), 5)
```

- `similar(v)` creates an uninitialized vector of the same type and size as `v`.

```{julia}
w = similar(v)
```



### Construction by Implicit Loop _(list comprehension)_

Implicit `for` loops are another method to create vectors.
```{julia}
v4 = [i for i in 1.0:8]
```


```{julia}
v5 = [log(i^2) for i in 1:4 ]
```

One can even insert an `if` clause.

```{julia}
v6 = [i^2 for i in 1:8 if i%3 != 2]
```

### Bit Vectors {#sec-bitvec}

Besides `Vector{Bool}`, there is also the special data type `BitVector` (and more generally `BitArray`) for storing arrays of truth values. 

While one byte is used for storing a `Bool`, storage in a BitVector is done bit by bit. 

The constructor converts a 
`Vector{Bool}` into a `BitVector`.

```{julia}
vb = BitVector([true, false, true, true])
```

For the reverse direction, there is `collect()`.

```{julia}
collect(vb)
```

Bit vectors are produced, for example, as the result of element-wise comparisons (s. @sec-broadcast).

```{julia}
v4 .> 3.5
``` 



### Indexing

Indices are ordinal numbers. Therefore, __indexing starts at 1.__

As an index, one can use:

 - Integer
 - Integer-valued range (same length or shorter)
 - Integer vector         (same length or shorter)
 - Boolean vector or BitVector (same length)
 
  
With indices, one can read and write array elements/parts.


```{julia}
v = [ 3i + 5.2 for i in 1:8]
```


```{julia}
v[5]
```

In assignments, the right side is converted to the vector element type with `convert(T,x)` if necessary.
```{julia}
v[6] = 9999
v
```


Exceeding the index limits leads to a `BoundsError`. 


```{julia}
v[77]
```

A `range` object can be used to address a subvector.

```{julia}
vp = v[3:5]
vp
```

```{julia}
vp = v[1:2:7]   # range with step size
vp
```

- When used as an index, the special value `end` can be used in a range.
- When used as an index, the "empty" range `:` can be used as an abbreviation for `1:end`. This is useful for matrices: `A[:, 3]` addresses the entire 3rd column of `A`. 


```{julia}
v[6:end] = [7, 7, 7]
v
```

#### Indirect Indexing

The indirect indexing with a *Vector of Integers/Indices* follows the formula  


$v[ [i_1,\ i_2,\ i_3,...]] = [\ v[i_1],\ v[i_2],\ v[i_3],...]$


```{julia}
v[ [1, 3, 4] ]
```

is therefore equal to
```{julia}
[ v[1], v[3], v[4] ]
```





#### Indexing with a Vector of Truth Values

As an index, one can also use a `Vector{Bool}` or `BitVector` (s. @sec-bitvec) **of the same length**. 

```{julia}
v[ [true, true, false, false, true, false, true, true] ]
```

This is useful as one can, for example,

 - broadcast tests (s. @sec-broadcast),
 - these tests then produce a BitVector and
 - Bit vectors can be combined with bitwise operators `&` and `|` as needed.    



```{julia}
v[  (v .> 13) .& (v.<20) ]
```




## Matrices and Arrays

The methods presented so far for vectors also apply to higher-dimensional arrays.

One can create them uninitialized:

```{julia}
A = Array{Float64,3}(undef, 6,9,3)
```

In most functions, the dimensions can also be passed as a tuple. The above instruction can also be written as:
```julia
A = Array{Float64, 3}(undef, (6,9,3))  
```

Functions like `zeros()` etc. of course also work.

```{julia}
m2 = zeros(3, 4, 2)  # or zeros((3,4,2))   
```


```{julia}
M = fill(5 , (3, 3))   # or fill(5, 3, 3)
```

The function `similar()`, which creates an uninitialized array of the same size, can also take a type as a further argument.

```{julia}
M2 = similar(M, Float64)
```


### Construction by Explicit Element List

While vectors are noted in square brackets separated by commas, the notation for higher-dimensional objects is somewhat different. 

- A matrix:

```{julia}
M2 = [2  3  -1
      4  5  -2]
```

- the same matrix:

```{julia}
M2 = [2 3 -1; 4 5 -2]
```

- An array with 3 indices:

```{julia}
M3 = [2   3  -1 
      4   5   6 ;;;
      7   8   9
     11  12  13]
```

- and once more the matrix `M2`:

```{julia}
M2 = [2;4;; 3;5;; -1;-2]
```

In the last example, these rules apply:

- The separator is the semicolon.
- A semicolon `;` increases the 1st index.
- Two semicolons `;;`  increase the 2nd index.
- Three semicolons `;;;` increase the 3rd index etc. 

In the previous examples, the following syntactic enhancement (_syntactic sugar_) was applied:

- Spaces separate like 2 semicolons -- thus also increase the 2nd index: $\quad a_{12}\quad a_{13}\quad a_{14}\ ...$
- Newline separates like a semicolon -- thus also increases the 1st index.


:::{.callout-important}

- Vector notation with comma as separator only works for vectors, do not mix "semicolon, space, newline"! 
- Vectors, $1\!\times\!n$-matrices, and $n\!\times\!1$-matrices are three different things!


```{julia}
v1 = [2,3,4]
```


```{julia}
v2 = [2;3;4]
```


```{julia}
v3 = [2 3 4]
```


```{julia}
v3 = [2;3;4;;]
```

:::


One can of course also construct a "vector of vectors" in the C/C++ style.  


```{julia}
v = [[2,3,4], [5,6,7,8]]
```


```{julia}
v[2][3]
```

One should only do this in special cases. The array language of Julia is usually more convenient and faster.

### Indices, Subarrays, Slices


```{julia}
# 6x6 matrix with random numbers uniformly distributed from [0,1) ∈ Float64
A = rand(6,6)
```

The usual index notation:

```{julia}
A[2, 3] = 77.77777
A
```

One can use ranges to address subarrays:

```{julia}
B = A[1:2, 1:3]
```

Addressing parts with lower dimension is also called *slicing*.

```{julia}
# the 3rd column as a vector (slicing)

C = A[:, 3]
```


```{julia}
# the 3rd row as a vector (slicing)

E = A[3, :]
```

Of course, assignments are also possible:

```{julia}
# One can also assign to slices and subarrays 

A[2, :] = [1,2,3,4,5,6]
A
```

## Behavior in Assignments, `copy()` and `deepcopy()`, Views

### Assignments and Copies

- Variables are references to objects. 
- _An assignment to a variable does not create a new object._



```{julia}
A = [1, 2, 3]
B = A
```

`A` and `B` are now names of the same object.

```{julia}
A[1] = 77
@show B;
```


```{julia}
B[3] = 300
@show A;
```

This behavior saves a lot of time and memory, but is not always desired. 
The function `copy()` creates a 'real' copy of the object.

```{julia}
A = [1, 2, 3]
B = copy(A)
A[1] = 100
@show A B;
```

The function 
`deepcopy(A)` copies recursively. Copies are also created for the elements from which `A` consists. 

As long as an array only contains primitive objects (numbers), `copy()` and `deepcopy()` are equivalent.  

The following example shows the difference between `copy()` and `deepcopy()`.
```{julia}
mutable struct Person
    name :: String
    age  :: Int
end

A = [Person("Meier", 20), Person("Müller", 21), Person("Schmidt", 23)]
B = A
C = copy(A)
D = deepcopy(A)
```

```{julia}
A[1] = Person("Mustermann", 83)
A[3].age = 199 

@show B C D;
```

### Views

When one assigns a piece of an array to a variable using *indices/ranges/slices*, 
Julia fundamentally **constructs a new object**.


```{julia}
A = [1  2  3
     3  4  5]

v = A[:, 2]
@show  v

A[1, 2] = 77
@show  A  v;
```


Sometimes, however, one wants exactly this reference semantics in the sense of: "Vector `v` should be the 2nd column vector of `A` and should also remain so (i.e., change if `A` changes)."

This is called *views* in Julia: We want the variable `v` to represent only an 'alternative view' of the matrix `A`.

This can be achieved with the `@view` macro: 

```{julia}
A = [1 2 3
     3 4 5]

v = @view A[:,2]
@show  v

A[1, 2] = 77
@show  v;
```


Julia uses this technique for efficiency reasons also in some functions of linear algebra. 
An example is the operator `'`, which delivers the adjoint matrix `A'` to a matrix `A`.

- The adjoint matrix `A'` is the transposed and element-wise complex-conjugated matrix to `A`. 
- The parser converts this to the function call `adjoint(A)`.  
- For real matrices, the adjoint is equal to the transposed matrix.
- Julia implements `adjoint()` as a _lazy function_, i.e., 
- for efficiency reasons, no new object is constructed, but only an alternative 'view' of the matrix ("with swapped indices") and an alternative 'view' of the entries (with sign change in the imaginary part).
- From vectors, `adjoint()` makes a $1\!\times\!n$-matrix (a row vector).


```{julia}
A = [ 1.  2.
      3.  4.]
B = A'
```

The matrix `B` is only a modified 'view' of `A`:

```{julia}
A[1, 2] =10
B
```


From vectors, `adjoint()` makes a $1\!\times\!n$-matrix (a row vector).


```{julia}
v = [1, 2, 3]
v'
```

Another such function, which provides an alternative 'view', a different indexing, of the same data 
is `reshape()`. 

Here, a vector with 12 entries is transformed into a 3x4 matrix.
```{julia}
A = [1,2,3,4,5,6,7,8,9,10,11,12]

B = reshape(A, 3, 4)   
```

## Storage of an Array

- Memory is addressed linearly. A matrix can be stored in memory row-wise _(row major)_ or column-wise _(column major)_. 
- C/C++/Python(NumPy) use row-major storage: The 4 elements of a 2x2 matrix are stored in the order $a_{11},a_{12},a_{21},a_{22}$. 
- Julia, Fortran, Matlab store column-wise: $a_{11},a_{21},a_{12},a_{22}$. 

This information is important to iterate efficiently over matrices:


```{julia}
function column_major_add(A, B)
    (n,m) = size(A)
    for j = 1:m
        for i = 1:n     # inner loop traverses a column
            A[i,j] += B[i,j]
        end
    end
end

function row_major_add(A, B)
    (n,m) = size(A)
    for i = 1:n
        for j = 1:m     # inner loop traverses a row
            A[i,j] += B[i,j]
        end
    end
end
```


```{julia}
A = rand(10000, 10000);
B = rand(10000, 10000);
```


```{julia}
using BenchmarkTools

@benchmark row_major_add($A, $B) 
```


```{julia}
@benchmark column_major_add($A, $B)
```

### Locality of Memory Access and _Caching_

We have seen that the order of inner and outer loops makes a significant speed difference:

__It is more efficient when the innermost loop traverses the left index__, i.e., a column and not a row. The cause of this lies in the architecture of modern processors.

- Memory access occurs over multiple cache levels. 
- A _cache miss_, which triggers a reload from slower caches, slows down the process.
- Larger memory blocks are always reloaded to minimize the frequency of _cache misses_.
- Therefore, it is important to organize memory access as locally as possible.

::: {.content-visible when-format="html"}
![Memory hierarchy of Intel processors, from: Victor Eijkhout,_Introduction to High-Performance Scientific Computing_, [https://theartofhpc.com/](https://theartofhpc.com/)](../images/cache.png){width="75%"}

:::

::: {.content-visible when-format="pdf"}
![Memory hierarchy of Intel processors, from: Victor Eijkhout,_Introduction to High-Performance Scientific Computing_, [https://theartofhpc.com/](https://theartofhpc.com/)](../images/cache.png){width="70%"}
:::


## Mathematical Operations with Arrays

Arrays of the same dimension (e.g., all $7\!\times\!3$-matrices) form a linear space. 

 - They can be multiplied by scalars and
 - they can be added and subtracted.


```{julia}
0.5 * [2, 3, 4, 5]
```


```{julia}
0.5 * [ 1  3
        2  7] - [ 2  3; 1 2]
```

### Matrix Product

The matrix product is defined for 

:::{.narrow}

| 1st factor | 2nd factor | Product |
| :-:  | :-:   | :-: |
| $(n\!\times\!m)$-matrix | $(m\!\times\!k)$-matrix | $(n\times k)$-matrix|
| $(n\!\times\!m)$-matrix | $m$-vector | $n$-vector |
| $(1\!\times\!m)$-row vector | $(m\!\times\!n)$-matrix | $n$-vector |
| $(1\!\times\!m)$-row vector | $m$-vector | scalar product |
| $m$-vector | $(1\times n)$-row vector | $(m\!\times\!n)$-matrix |

:::

Examples:
```{julia}
A = [1 2 3
     4 5 6]
v = [2, 3]
w = [1, 3, 4];
```

- (2,3)-matrix `*` (3,2)-matrix
```{julia}
A * A'
```

- (3,2)-matrix `*` (2,3)-matrix

```{julia}
A' * A
```

- (2,3)-matrix `*` 3-vector

```{julia}
A * w
```

- (1,2)-vector `*` (2,3)-matrix

```{julia}
v' * A
```

- (3,2)-matrix `*` 2-vector

```{julia}
A' * v
```

- (1,2)-vector `*` 2-vector (scalar product)

```{julia}
v' * v
```

2-vector `*` (1,3)-vector (outer product) 

```{julia}
v * w'
```

## Broadcasting {#sec-broadcast}

- With _broadcasting_, operations or functions are applied __element-wise__ to arrays.
- The syntax for this is a dot _before_ an operator or _after_ a function name. 
- The parser converts `f.(x,y)` to `broadcast(f, x, y)` and similarly for operators `x .⊙ y` to `broadcast(⊙, z, y)`. 
- Operands that are missing one or more dimensions are virtually replicated in those dimensions. 
- The *broadcasting* of assignments `.=`, `.+=`,... changes the semantics. No new object is created, but the values are entered into the object on the left side (which must have the correct dimension).

 

Some examples:

- Element-wise application of a function
```{julia}
sin.([1, 2, 3])
```

- The following does not give the algebraic [square root of a matrix](https://en.wikipedia.org/wiki/Square_root_of_a_matrix), but the element-wise square root of each entry. 
```{julia}
A = [8 2
     3 4]

sqrt.(A)
```

- The following does not give $A^2$, but the entries are squared.
```{julia}
A.^2   
```

- For comparison, the result of the algebraic operations:
```{julia}
@show A^2 A^(1/2);
```

- Broadcasting also works with functions of several variables.
```{julia}
hyp(a,b) = sqrt(a^2+b^2)

B = [3 4
     5 7]

hyp.(A, B) 
```



When operands have different dimensions, the operand with missing dimensions is virtually 
'grown' by replication in these dimensions.

We add a scalar to a matrix:
```{julia}
A = [ 1 2 3
      4 5 6]
```


```{julia}
A .+ 300
```

The scalar was replicated to the same dimension as the matrix by replication. We let `broadcast()` show the form of the 2nd operand after *broadcasting*: 

```{julia}
broadcast( (x,y) -> y, A, 300)
```

(Of course, this replication only takes place virtually. This object is not really created for other operations.)

As another example: Matrix and (column-)vector

```{julia}
A .+ [10, 20]
```

The vector is expanded by repeating columns:

```{julia}
broadcast((x,y)->y, A, [10,20])
```

Matrix and row vector: The row vector is repeated row-wise:

```{julia}
A .* [1,2,3]'      # adjoint vector
```

The 2nd operand is 'grown' by `broadcast()` through replication of rows.

```{julia}
broadcast((x,y)->y, A, [1,2,3]')
```


#### _Broadcasting_  in Assignments

Assignments  `=`, `+=`, `/=`,..., where a name is on the left side, work in Julia such that an object is constructed from the right side and assigned this new name. 

However, when working with arrays, one often wants to reuse an existing array for efficiency reasons. The values calculated on the right side should be entered into the already existing object on the left side. 

This is achieved with the broadcast variants  `.=`, `.+=`,... of the assignment operators.

```{julia}
A .= 3
```

```{julia}
A .+= [1, 4]
```


## Further Array Functions - a Selection

Julia provides a large number of functions that work with arrays.

```{julia}
A = [22 -17 8 ; 4 6 9]
```

- Find the maximum

```{julia}
maximum(A)
```

- Find the maximum of each column
```{julia}
maximum(A, dims=1)  
```

- Find the maximum of each row
```{julia}
maximum(A, dims=2)
```

- Find the minimum and its position
```{julia}
amin, i = findmin(A)
```

- What is a `CartesianIndex`?
```{julia}
dump(i)
```

- Extract the indices of the minimum as a tuple
```{julia}
i.I
```

- Sum and product of all entries
```{julia}
sum(A), prod(A)
```

- Column sum (1st index is reduced)

```{julia}
sum(A, dims=1) 
```

 
- Row sums   (2nd index is reduced)

```{julia}
sum(A, dims=2) 
```

- Sum after applying a function element-wise

```{julia}
sum(x->sqrt(abs(x)), A)   #   sum_ij sqrt(|a_ij|)
```

- Reduce (fold) the array with a function

```{julia}
reduce(+, A)    #  equivalent to sum(A)
```

-  `mapreduce(f, op, array)`: Apply `f` to all entries, then reduce with `op`

```{julia}

mapreduce(x -> x^2, +, A )     # Sum of the squares of all entries
```


- Are there elements in A that are > 5? 

```{julia}
any(x -> x>5, A)    
```

- How many elements in A are > 5?

```{julia}
count(x-> x>5, A)  
```

-  are all entries positive?
```{julia}
all(x-> x>0, A)
```

 
