---
engine: julia
---

```{julia}
#| error: false
#| echo: false
#| output: false
using InteractiveUtils

#struct M a::Int end; x = M(22); @show x 
#should not print "Main.Notebook.M(22)" but only "M(22)"
function Base.show(io::IO, x::T) where T
    if parentmodule(T) == @__MODULE__
        # Print "TypeName(fields...)" without module prefix
        print(io, nameof(T), "(")
        fields = fieldnames(T)
        for (i, f) in enumerate(fields)
            print(io, getfield(x, f))
            i < length(fields) && print(io, ", ")
        end
        print(io, ")")
    else
        invoke(Base.show, Tuple{IO, Any}, io, x)
    end
end
```

# Vectors, Matrices, Arrays

## General 

We now turn to containers important for numerical computing: 

- 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 from an explicit element list, Julia determines the greatest common type for the type parameter `T`.

```{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 'classical' arrays for scientific computing is that they are simply contiguous memory segments where same-length components (e.g., 64 bits) are stored sequentially. This minimizes memory requirements and maximizes access speed for both reading and modifying components. The location of `v[i]` can be calculated directly from `i`. 
- Julia's `Array{T,N}` (including vectors and matrices) is implemented this way for standard numeric types `T`. Elements are stored *unboxed*. In contrast, `Vector{Any}` is implemented as a list of object addresses *(boxed)*, 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

### Stack-like Functions 

- `push!(vector, items...)`   --- appends elements to the end
- `pushfirst!(vector, items...)` --- prepends elements to the beginning
- `pop!(vector)` --- removes and returns the last element
- `popfirst!(vector)` --- removes and returns the first element  

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

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

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

push!(v, 17)
```

A `push!()` operation can be expensive, as it may require allocating new memory and copying the existing vector. Julia optimizes memory by preallocating space, so subsequent `push!` operations are very fast, almost achieving O(1) speed.   

Avoid `push!()` or `resize()` in time-critical code with very large arrays.

### Further Constructors

You can create uninitialized vectors of a given length and type. This is the fastest method; elements contain random bit patterns.
```{julia}
# Uninitialized vector of length 1000

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


- `zeros(n)` creates a `Vector{Float64}` of length `n`, initialized with zeros.  

```{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`, filled 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 provide another method to create vectors.
```{julia}
v4 = [i for i in 1.0:8]
```


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

You can even include an `if` clause.

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

### Bit Vectors {#sec-bitvec}

Besides `Vector{Bool}`, Julia provides the `BitVector` data type (and more generally `BitArray`) for storing boolean arrays. 

While a `Bool` uses one byte, `BitVector` stores values bit by bit. 

The constructor converts a `Vector{Bool}` to a `BitVector`:

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

Use `collect()` for the reverse conversion:

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

Bit vectors are produced, for example, by element-wise comparisons (see @sec-broadcast):

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



### Indexing

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

Valid indices include:

 - an integer
 - an integer-valued range (same length or shorter)
 - an integer vector         (same length or shorter)
 - a boolean vector or BitVector (same length)
 
  
Using indices, we 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 index limits throws 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

Indirect indexing with a *vector of integers* follows the formula


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


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

which is the same as
```{julia}
[ v[1], v[3], v[4] ]
```





#### Indexing with a Vector of Truth Values

You can also use a `Vector{Bool}` or `BitVector` (see @sec-bitvec) **of the same length** as an index. 

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

This is useful for:

 - broadcast tests (see @sec-broadcast) which produce a `BitVector`, and
 - combining `BitVector`s with bitwise operators `&` and `|`.    


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




## Matrices and Arrays

Most methods 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 can be written as:
```julia
A = Array{Float64, 3}(undef, (6,9,3))  
```

Functions like `zeros()` etc. 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 written 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]
```

Here, the following 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 two semicolons -- thus increasing the 2nd index: $\quad a_{12}\quad a_{13}\quad a_{14}\ ...$
- Newline separates like a semicolon -- thus increasing 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\times1$-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]
```

You should only do this in special cases. The array notation in 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
```

You 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, :]
```

Slicing can also used on the right hand side for  assignments:

```{julia}
# You 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 true copy of the object.

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

The function 
`deepcopy(A)` copies recursively. It also creates copies of the elements that A contains. 

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  **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  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 a *view* in Julia: We want the variable `v` to represent only an 'alternative view' of the matrix `A`.

It 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. The method provides an alternative 'view' of the matrix (with swapped indices) and an alternative 'view' of the entries (with sign change in the imaginary part).
- The adjoint of a vector produces a $1\times n$ matrix (row vector).


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

The matrix `B` is just 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 observed that the order of inner and outer loops significantly affects computational efficiency:

It is more efficient when the innermost loop iterates over the left index, i.e., a column rather than a row. This is due to the architecture of modern processors.

- Memory access involves multiple cache levels.
- A _cache miss_, which necessitates reloading from slower caches, slows down processing.
- To minimize the frequency of _cache misses_, processors reload larger memory blocks.
- Consequently, it is crucial 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="typst"}
![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\times3$-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.
- Broadcasting is indicated by a dot preceding an operator or following a function name.
- The parser translates `f.(x,y)` into `broadcast(f, x, y)`, and similarly, `x .⊙ y` into `broadcast(⊙, x, y)`.
- Operands lacking one or more dimensions are virtually replicated in those dimensions.
- Broadcasting assignments such as `.=`, `.+=`, etc., alter the semantics by modifying values directly within the left-side object (which must have the appropriate dimensions) without creating a new object.
- With _broadcasting_, operations or functions are applied element-wise to arrays.

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, here the results 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 possess differing dimensions, the operand with fewer dimensions is effectively expanded through replication along those dimensions.

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


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

The scalar was replicated to match the matrix dimensions. Let `broadcast()` illustrate the form of the second operand after broadcasting:
```
broadcast( (x,y) -> y, A, 300)
```

(This replication occurs solely in a virtual sense; the object is not actually instantiated in memory.)

```{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 such as `=`, `+=`, `/=`,... in Julia involve constructing an object on the right-hand side and assigning it to a variable on the left-hand side.

When working with arrays, efficiency often requires reusing existing array objects. The values computed on the right-hand side are then stored directly into the pre-existing object on the left-hand side.

This is accomplished using broadcast variants of 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)
```

 
