attempt to explain why broadcasting is bae

broadcast.Rmd

A starter example is centering a matrix

```{python}
import numpy as np

X = np.random.normal(size=(3, 3))  # random 3 by 3 matrix

mu = X.mean(axis = 0)    # array (vector) of column means

print("Original matrix")
print(X)

print("\nArray of column means")
print(mu)

# array effectively copied till vector matches matrix in size
print("\nBroadcast subtraction operation")
print(X - mu)  
```

To do this in R you have to loop or use something like map/apply.

This lets you create kernel matrices in a fast and fairly readable way, for example.

```{python}
import numpy as np

X = np.random.normal(size=(10, 5))
Y = np.random.normal(size=(2, 5))
sigma = 0.1

# calculate all pairwise distances between rows in X and Y 
# (x - y)^2 = a^2 - 2 ab + b^2
dist_sq = np.square(X).sum(axis=1)[:, np.newaxis] - 2.0 * np.dot(X, Y.T) + np.square(Y).sum(axis=1)

# turn this into an RBF kernel (math might be off, but hopefully the point is clear)
rbf = np.exp(-np.sqrt(dist_sq) / (2 * sigma ** 2))

print(dist_sq)
print(rbf)
```

In R you have to hope someone vectorized the components of the operation that
you care about, copy an object along an appropriate dimension, or use
loops/applys/maps, which I think make the math harder to think about. 

Just discovered this, don't know how I hadn't seen your comment yet -- super cool!

FYI if you dig broadcasting, the whole point of xtensor is to provide the same Numpy interface and broadcasting niceness from within Rcpp.

Ugh if only that played better with windows.

So this is how rray came to be hmm 🤔