Hi Bradley,

Sorry, my explanation of that hack was brief, mostly because it was
only a hack.  The key comment was:

  "However, I intend to flatten all of the affine transformations
next, at which point I'll need to implement this properly."

The entire point of composing affine transformations is that their
composition is another affine transformation, so the indexer never
grows in complexity.  In particular, for a 2D matrix this
transformation is just a matter of strides, which means BLAS can
multiply any SRFI 179 f32 or f64 vectors.

What surprised me most is that after completing my work and looking at
the reference implementation, the latter does not seem to perform this
optimization.  In fact, the array record type only holds the indexer,
not the affine coefficients, so it cannot perform this optimization
(unless the compiler is smart enough to combine these closures at
runtime).

--
Alex

On Sat, Apr 17, 2021 at 8:46 AM Bradley Lucier <xxxxxx@math.purdue.edu> wrote:
>
> Dear Alex:
>
> Your routine
>
> (define (specialized-array-share array new-domain project)
>    (assert (specialized-array? array) (interval? new-domain))
>    (let ((body (array-body array))
>          (indexer (lambda multi-index
>                     (call-with-values
>                         (lambda () (apply project multi-index))
>                       (array-indexer array))))
>          (storage (array-storage-class array)))
>      (%make-array
>       new-domain
>       (specialized-getter body indexer (storage-class-getter storage))
>       (specialized-setter body indexer (storage-class-setter storage))
>       storage
>       body
>       indexer
>       (array-safe? array))))
>
> does not take into account the affine structure of either project or
> (array-indexer array) to simplify the indexer of the result array as in
> the procedure (array/optimize-linear-function f d) in Bawden's original
> newsgroup post
>
> https://groups.google.com/g/comp.lang.scheme/c/7nkx58Kv6RI/m/a5hdsduFL2wJ
>
> I thought to look more carefully after I installed the git version of
> chibi myself to try the example in my previous email and everything worked:
>
> =========================================================
> heine:~/programs/chibi-scheme> /usr/local/chibi/bin/chibi-scheme
>  > (import (srfi 179))
>  > (define A (array-copy (make-array (make-interval '#(3 4)) list)))
>  > (define B (array-sample A '#(2 1))) >
>  > (define (array-display A)
>
>    (define (display-item x)
>      (display x) (display "\t"))
>
>    (newline)
>    (case (array-dimension A)
>      ((1) (array-for-each display-item A) (newline))
>      ((2) (array-for-each (lambda (row)
>                             (array-for-each display-item row)
>                             (newline))
>                           (array-curry A 1)))
>      (else
>       (error "array-display can't handle > 2 dimensions: " A))))
>  > (array-display (specialized-array-reshape B (make-interval '#(8))))
>
> (0 0)   (0 1)   (0 2)   (0 3)   (0 0)   (0 1)   (0 2)   (0 3)
>  > (array-sample (specialized-array-reshape B (make-interval '#(8))) '#(3))
> #<Array 139698477832576>
>  > (array-display (array-sample (specialized-array-reshape B
> (make-interval '#(8))) '#(3)))
>
> (0 0)   (0 3)   (0 2)
>  > (array-body B)
> #((0 0) (0 1) (0 2) (0 3) (1 0) (1 1) (1 2) (1 3) (2 0) (2 1) (2 2) (2 3))
>  > (array-display A)
>
> (0 0)   (0 1)   (0 2)   (0 3)
> (1 0)   (1 1)   (1 2)   (1 3)
> (2 0)   (2 1)   (2 2)   (2 3)
>  > (array-display B)
>
> (0 0)   (0 1)   (0 2)   (0 3)
> (2 0)   (2 1)   (2 2)   (2 3)
>  > (specialized-array? (array-sample (specialized-array-reshape B
> (make-interval '#(8))) '#(3)))
> #t
> =========================================================
>
> I think implicit in the text of this SRFI is that the indexer of a
> specialized array should take time bounded independent of the number of
> Bawden-style transformations that have been applied to it.  In the SRFI
> document there are the following two statements:
> =========================================================
> Second, because the composition of any number of affine mappings is
> again affine, accessing or changing the elements of a restricted,
> translated, curried, permuted array is no slower than accessing or
> changing the elements of the original array itself.
> =========================================================
> where those operations are applied in order to the same array, and
> =========================================================
> Note: It is assumed that the affine structure of the composition of
> new-domain->old-domain and (array-indexer array) will be used to simplify:
>
>    (lambda multi-index
>      (call-with-values
>          (lambda ()
>            (apply new-domain->old-domain multi-index))
>        (array-indexer array)))
> =========================================================
>
> I think using your definition of specialized-array-reshape in the sample
> implementation, while leaving everything else the same, would cause
> =========================================================
>  > (array-display (array-sample (specialized-array-reshape B
> (make-interval '#(8))) '#(3)))
>
> (0 0)   (0 3)   (0 2)   =========================================================
> from the above example to fail.
>
> I don't know whether more explicit language is needed in the SRFI
> document to guide implementors.
>
> Brad