I'm trying to understand array broadcasting and for what purposes is it
used.  I'm soliciting the perspectives of others.

TL;DR summary: If the procedures in SRFI 231 offer substantially the
same power as array broadcasting, perhaps we don't need to have array
broadcasting in Scheme.

===========================

I searched the web for example code in NumPy, but because array
broadcasting is implicitly invoked, it's a bit hard to Google for "real
world" examples, as it were.

So I worked through some examples in tutorials and old NumPy
documentation I found on the web using what's currently available in
SRFI 231.

Because of this exercise, it appears to me that there are (at least) two
different styles of array programming in general.

One is APL style, which SRFI 231 imitates in some ways with array-map,
array-curry, array-decurry, array-permute, array-outer-product,
array-inner-product, ...

NumPy doesn't particularly use this style of programming.  Instead, it
has array broadcasting (together with inserting new axes) as a sort of
"one weird trick" that allows to imitate the APL-type procedures.

If people are used to to thinking in this NumPy style of programming, or
if Scheme programmers think NumPy's style of array programming is
important, then we should try to implement some version of it.

But if the procedures in SRFI 231 offer substantially the same power as
array broadcasting, perhaps we don't need to have array broadcasting in
Scheme.

Worked examples follow.

Brad

(declare
   (standard-bindings)
   (extended-bindings)
   (block)
   (not safe))

(import (srfi 231))

;;; Examples from
;;;
https://www.pythonlikeyoumeanit.com/Module3_IntroducingNumpy/Broadcasting.html

;;; Utilizing Size-1 Dimensions for Broadcasting

(define images
   (let ((domain
          (make-interval
           '#(500   ;; number of images
              48    ;; number of rows/image
              48    ;; number of columns/image
              3)))) ;; three channels per pixel
     (specialized-array-reshape
      (make-specialized-array-from-data
       (random-f64vector (interval-volume domain))
       f64-storage-class)
      domain)))

(define channel-max
   (array-copy
    (array-map
     (lambda (channel)
       (array-reduce max channel))
     (array-curry (array-permute images '#(0 3 1 2)) 2))))

(define scaled-images
   (array-permute
    (array-decurry
     (array-map (lambda (channel channel-max)
                  (array-map (lambda (pixel)
                               (/ pixel channel-max))
                             channel))
                (array-curry (array-permute images '#(0 3 1 2)) 2)
                channel-max))
    '#(0 2 3 1)))

;;; An Advanced Application of Broadcasting: Pairwise Distances

(define (make-random-3d-points number-of-points)
   (let ((domain (make-interval (vector number-of-points 3))))
     (specialized-array-reshape
      (make-specialized-array-from-data
       (random-f64vector (interval-volume domain))
       f64-storage-class)
      domain)))

(define x-points (time (make-random-3d-points 10000)))
(define y-points (time (make-random-3d-points 10000)))

(define (Frobenius A B)
   (flsqrt (array-reduce fl+ (array-map (lambda (a b) (flsquare (fl- a
b))) A B))))

(define mutual-distances
   (time
    (array-copy!
     (array-outer-product
      Frobenius
      (array-copy! (array-curry x-points 1))
      (array-copy! (array-curry y-points 1))))))

(array->list* (array-extract mutual-distances (make-interval '#(20 5))))

;;; From
;;; https://scipy.github.io/old-wiki/pages/EricsBroadcastingDoc

;;; Comment at the end of this page:

#|
Broadcasting is a powerful tool for writing short and usually intuitive
code that does its computations very efficiently in C. However, there
are cases when broadcasting uses unnecessarily large amounts of memory
for a particular algorithm. In these cases, it is better to write the
algorithm’s outer loop in Python. This may also produce more readable
code, as algorithms that use broadcasting tend to become more difficult
to interpret as the number of dimensions in the broadcast increases.
|#

;;; Vector quantization example (references at page)

(define number-of-observations 4000)

(define number-of-codes 40)

(define number-of-features 16)

(define observations
   (let ((domain (make-interval (vector number-of-observations
number-of-features))))
     (specialized-array-reshape
      (make-specialized-array-from-data
       (random-f64vector (interval-volume domain))
       f64-storage-class)
      domain)))

(define codes
   (let ((domain (make-interval (vector number-of-codes
number-of-features))))
     (specialized-array-reshape
      (make-specialized-array-from-data
       (random-f64vector (interval-volume domain))
       f64-storage-class)
      domain)))

(define curried-codes
   (array-copy (array-curry codes 1)))

(define (argmin f interval)
   (cdr (interval-fold-left (lambda args
                              (cons (apply f args) args))
                            (lambda (existing new)
                              (if (< (car new)
                                     (car existing))
                                  new existing))
                            '(+inf.0) ;;; don't need any args
                            interval)))

(define (Frobenius A B)
   (flsqrt (array-reduce fl+ (array-map (lambda (a b) (flsquare (fl- a
b))) A B))))

(define code-vectors
   (time
    (array-copy (array-map (lambda (observation)
                             (car  ;; the result is a list of length 1,
we just need the index
                              (argmin (lambda (i)
                                        (Frobenius observation
(array-ref curried-codes i)))
                                      (array-domain curried-codes))))
                           (array-curry observations 1)))))

(define (bin array vector)
   (array-for-each (lambda (i)
                     (vector-set! vector i (+ 1 (vector-ref vector i))))
                   array)
   vector)

#|
(bin code-vectors (make-vector number-of-codes 0))
#(36 126 193 351 82 145 71 59 143 124 74 92 52 92 81 113 39 54 124 20 19
104 101 43 152 77 69 249 112 74 25 111 65 38 146 128 158 67 103 88)
|#

#|
(time (array-copy (array-map (lambda (observation) (car (argmin (lambda
(i) (Frobenius observation (array-ref curried-codes i))) (array-domain
curried-codes)))) (array-curry observations 1))))
     0.151801 secs real time
     0.151092 secs cpu time (0.151092 user, 0.000000 system)
     no collections
     215586096 bytes allocated
     no minor faults
     no major faults
     545186244 cpu cycles
|#

;;; OpenXLA has fewer implicit conversions when broadcasting:

;;; https://openxla.org/xla/broadcasting