---
title: "Advanced indexing & Broadcasting"
subtitle: "Lecture 06"
author: "Dr. Colin Rundel"
footer: "Sta 663 - Spring 2023"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute:
  echo: true
---

```{r, message=FALSE, warning=FALSE, include=FALSE}
options(
  width=80
)

local({
  hook_old <- knitr::knit_hooks$get("error")  # save the old hook
  knitr::knit_hooks$set(error = function(x, options) {
    x = sub("## \n## Detailed traceback:\n.*$", "", x)
    x = sub("Error in py_call_impl\\(.*?\\)\\: ", "", x)
    hook_old(x, options)
  })
})
```

```{python include=FALSE}
import numpy as np

np.set_printoptions(edgeitems=3, linewidth=180)
```



# Advanced Indexing

## Advanced Indexing

> Advanced indexing is triggered when the selection object, obj, is a non-tuple sequence object, an ndarray (of data type integer or bool), or a tuple with at least one sequence object or ndarray (of data type integer or bool).
> 
> * There are two types of advanced indexing: integer and Boolean.
>
> * Advanced indexing always returns a *copy* of the data (contrast with basic slicing that returns a view).


::: {.aside}
From NumPy - User Guide - [Advanced Indexing](https://numpy.org/doc/stable/user/basics.indexing.html#advanced-indexing)
::: 

## From last time: subsetting with tuples

Unlike lists, an ndarray can be subset by a tuple containing integers,

:::: {.columns .small}
::: {.column width='50%'}
```{python}
x = np.arange(16).reshape((4,4)); x
x[(0,1,3), :]
x[:, (0,1,3)]
np.shares_memory(x, x[(0,1,3), :])
```
:::

::: {.column width='50%' .fragment}
```{python error=TRUE}
x[(0,1,3),]
x[(0,1,3)]
x[(0,1)]
x[0,1]
```
:::
::::



## Integer array subsetting (lists)

Lists of integers can be used to subset in the same way:

::: {.small}
```{python error=TRUE}
x = np.arange(16).reshape((4,4)); x
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python error=TRUE}
x[[1,3]]
x[[1,3], ]
x[:, [1,3]]
```
:::

::: {.column width='50%'}
```{python error=TRUE}
x[[0,1,3],]
x[[0,1,3]]
x[[1.,3]]
```
:::
::::


## Integer array subsetting (ndarrays)

Similarly we can also us integer ndarrays:

:::: {.columns .small}
::: {.column width='50%'}
```{python error=TRUE}
x = np.arange(6)
y = np.array([0,1,3])
z = np.array([1., 3.])
```

```{python error=TRUE}
x[y,]
x[y]

x[z]
```
:::

::: {.column width='50%' .fragment}
```{python}
x = np.arange(16).reshape((4,4))
y = np.array([1,3])
```

```{python}
x[y]
x[y, ]
x[:, y]
x[y, y]
```
:::
::::


## Exercise 1

Given the following matrix,

```{python}
x = np.arange(16).reshape((4,4))
x
```

write an expression to obtain the center 2x2 values (i.e. 5, 6, 9, 10 as a matrix).

## Boolean indexing

Lists or ndarrays of boolean values can also be used to subset (positions with `True` are kept and `False` are discarded)

::: {.small}
```{python}
x = np.arange(6); x
```
:::

. . .

::: {.small}
```{python error=TRUE}
x[[True, False, True, False, True, False]]
x[[True]]
```
:::

. . .

::: {.small}
```{python error=TRUE}
x[np.array([True, True, False, False, True, False])]
x[np.array([True])]
```
:::


## Boolean  expressions

The primary utility of boolean subsetting comes from vectorized comparison operations,

::: {.small}
```{python}
x = np.arange(6); x
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
x > 3
x[x>3]

x % 2 == 1
x[x % 2 == 1]
```
:::

::: {.column width='50%' .fragment}
```{python}
y = np.arange(9).reshape((3,3))
y % 2 == 0
y[y % 2 == 0]
```
:::
::::


## NumPy and Boolean operators

If we want to use a logical operators on an array we need to use `&`, `|`, and `~` instead of `and`, `or`, and `not` respectively.

::: {.medium}
```{python}
x = np.arange(6); x

y = x % 2 == 0; y

~y

y & (x > 3)

y | (x > 3)
```
:::


## meshgrid()

One other useful function in NumPy is `meshgrid()` which generates all possible combinations between the input vectors (as a tuple of ndarrays),

::: {.medium}
```{python}
pts = np.arange(3)
x, y = np.meshgrid(pts, pts)
x
y

np.sqrt(x**2 + y**2)
```
:::

## Exercise 2

We will now use this to attempt a simple brute force approach to numerical optimization, define a grid of points using `meshgrid()` to approximate the minima the following function:

$$
f(x,y) = (1-x)^2 + 100(y-x^2)^2
$$ 
Considering values of $x,y \in (-1,3)$, which value(s) of $x,y$ minimize this function?


```{r echo=FALSE, message=FALSE, out.width="60%", fig.align="center"}
library(tidyverse)
expand.grid(
    x = seq(-1,3, length.out = 1000),
    y = seq(-1,3, length.out = 1000)
) %>%
    mutate(
        f = (1-x)^2 + 100*(y-x^2)^2
    ) %>%
    ggplot(aes(x=x,y=y,fill=log(f))) +
        geom_raster() +
        theme_minimal() +
        labs(title = "Rosenbrock function (a=1, b=100)") +
        coord_fixed()
```



# Broadcasting

## Broadcasting

This is an approach for deciding how to generalize arithmetic operations between arrays with differing shapes.

```{python}
x = np.array([1, 2, 3])
x * 2
x * np.array([2])
x * np.array([2,2,2])
```

In the first example `2` is equivalent to the array `np.array([2])` which is being broadcast across the longer array `x`.

## Efficiancy

Using broadcasts can be much more efficient as it does not copy the underlying data,

```{python eval = FALSE}
x = np.arange(1e5)
y = np.array([2]).repeat(1e5)
```

```python
%timeit x * 2
17.3 µs ± 101 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
```

```python
%timeit x * np.array([2])
17.2 µs ± 239 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
```

```python
%timeit x * y
36 µs ± 121 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
```

## General Broadcasting

> When operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing (i.e. rightmost) dimensions and works its way left. Two dimensions are compatible when
>
>    1. they are equal, or
>
>    2. one of them is 1
> 
> If these conditions are not met, a `ValueError: operands could not be broadcast together` exception is thrown, indicating that the arrays have incompatible shapes. The size of the resulting array is the size that is not 1 along each axis of the inputs.

## Example

Why does the code on the left work but not the code on the right?

:::: {.columns .small}
::: {.column width='50%'}
```{python}
x = np.arange(12).reshape((4,3))
x

x + np.array([1,2,3])
```
:::

::: {.column width='50%'}
```{python error=TRUE}
x = np.arange(12).reshape((3,4))
x

x + np.array([1,2,3])
```
:::
::::

. . .

:::: {.columns}
::: {.column width='50%'}
```
x    (2d array): 4 x 3
y    (1d array):     3 
----------------------
x+y  (2d array): 4 x 3
```
:::

::: {.column width='50%'}
```
x    (2d array): 3 x 4
y    (1d array):     3 
----------------------
x+y  (2d array): Error
```
:::
::::



## A quick fix

```{python error=TRUE}
x = np.arange(12).reshape((3,4)); x
x + np.array([1,2,3]).reshape(3,1)
```

. . .

```
x    (2d array): 3 x 4
y    (2d array): 3 x 1
----------------------
x+y  (2d array): 3 x 4
```

## Examples (2)

:::: {.columns}
::: {.column width='50%'}
```{python}
x = np.arange(12).reshape((4,3))
y = 1
x+y
```

```
x    (2d array): 4 x 3
y    (1d array):     1 
----------------------
x+y  (2d array): 4 x 3
```
:::

::: {.column width='50%' .fragment}
```{python}
x = np.arange(12).reshape((4,3))
y = np.array([1,2,3])
x+y
```


```
x    (2d array): 4 x 3
y    (1d array):     3 
----------------------
x+y  (2d array): 4 x 3
```
:::
::::


## Examples (3)

::: {.small}
```{python}
x = np.array([0,10,20,30]).reshape((4,1))
y = np.array([1,2,3])
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
x
y
```
:::

::: {.column width='50%'}

```{python}
x+y
```
:::
::::

. . .

![](imgs/numpy_broadcasting.png){fig-align="center" width="50%"}

::: {.aside}
From NumPy user guide - [Broadcasting](https://numpy.org/doc/stable/user/basics.broadcasting.html)
:::

## Exercise 3

For each of the following combinations determine what the resulting dimension will be using broadcasting

* A [128 x 128 x 3] + B [3]

* A [8 x 1 x 6 x 1] + B [7 x 1 x 5]

* A [2 x 1] + B [8 x 4 x 3]

* A [3 x 1] + B [15 x 3 x 5]

* A [3] + B [4]


## Demo 1 - Standardization

Below we generate a data set with 3 columns of random normal values. Each column has a different mean and standard deviation which we can check with `mean()` and `std()`.

```{python}
rng = np.random.default_rng(1234)
d = rng.normal(loc=[-1,0,1], scale=[1,2,3], size=(1000,3))
```

```{python}
d.mean(axis=0)
d.std(axis=0)
```

Use broadcasting to standardize all three columns to have mean 0 and standard deviation 1. 

Check the new data set using `mean()` and `std()`.


## Broadcasting and assignment

In addition to arithmetic operators, broadcasting can be used with assignment via array indexing,

```{python error=TRUE}
x = np.arange(12).reshape((3,4))
y = -np.arange(4)
z = -np.arange(3)
```

:::: {.columns}
::: {.column width='50%'}
```{python}
x[:] = y
x
x[...] = y
x
```
:::

::: {.column width='50%'}
```{python error=TRUE}
x[:] = z
x[:] = z.reshape((3,1))
x
```
:::
::::



# Basic file IO

## Reading and writing ndarrays

We will not spend much time on this as most data you will encounter is more likely to be a tabular format (e.g. data frame) and tools like Pandas are more appropriate.

For basic saving and loading of NumPy arrays there are the `save()` and `load()` functions which use a built in binary format.

```{python}
x = np.arange(1e5)

np.save("data/x.npy", x)

new_x = np.load("data/x.npy")

np.all(x == new_x)
```

Additional functions for saving (`savez()`, `savez_compressed()`, `savetxt()`) exist for saving multiple arrays or saving a text representation of an array.

## Reading delimited data

While not particularly recommended, if you need to read delimited (csv, tsv, etc.) data into a NumPy array you can use `genfromtxt()`,

```{r include=FALSE}
options(width=300)
```

::: {.small}
```{python}
with open("data/mtcars.csv") as file:
    mtcars = np.genfromtxt(file, delimiter=",", skip_header=True)
    
mtcars
```
:::