dyadic transpose, a personal history

[Roger|Dyalog]

Over the years I have been thinking about ⍉ and have written about it.

In the following, the transpose functions are renamed tr79, tr16a, tr17, etc. from the originals in order to facilitate reference and comparisons.

1977-1981

The July 1977 issue of the Hewlett-Packard Journal [HP 1977] dealt exclusively with APL. In that issue, the article APL Data: Virtual Workspaces and Shared Storage by Grant Munsey describes the subscript calculus in Phil Abrams’s Ph.D. thesis [Abrams 1970], which relates row-major access of an array (equation 1, p.7)



to an alternative access (equation 2, p. 8)



Munsey wrote that the function take, drop, reversal, transpose, and reshape can be implemented in this alternative access without moving or copying the data.

At the time, I did not have access to the Abrams thesis. The relevant section is on page 73:

      d. a⍉m

      r ← rvec
      d ← del
      rank ← 1+(⌈/a)
      i ← 0
      del ← rank↑del
      rvec ← rank↑rvec
      rank   repeat
        begin
          rvec[i] ← ⌊/(i=a)/r
          del[i] ← +/(i=a)/d
          i ← i+1
      end

It was not difficult to devise a model of x⍉y from equation 2 even without access to the Abrams thesis. ⍴x⍉y was simple enough and was my one contribution to Paul Berry’s SHARP APL Reference Manual, on page 146:

      (⍴z)   ←→ (⍴y)⌊.+(⌈/⍴y)×x∘.≠⍳(⌈/x)+~⎕io    ⍝ 1979-03
      (⍴x⍉y) ←→ (⍴y)⌊.+(⌈/⍴y)×x∘.≠⍳⌈/0,x+~⎕io    ⍝ 1981-06

I don’t have a record of the 1-line model. It is an interesting exercise in “archæological algorithmic analysis” to reconstruct what it could have been in 1979. From the 2016 section below, we have

      tr16a ← {(,⍵)⌷⍨⊂(↑,¨⍳(⍴⍵)⌊.+(⌈/⍴⍵)×~b)+.×(⍴⍵)⊥b←⍺∘.=⍳0⌈1+⌈/⍺}

itself a forward reconstruction of the 1987 version (1987 section below) into current Dyalog APL. Whence:

      z←x tr79 y;l;s;t                                     
      z←(,y)[(s⊥l)+.×t⊤t⍴⍳×/t←s⌊.+(⌈/s←⍴y)×~l←x∘.=⍳0⌈1+⌈/x]

where:
  t←s⌊.+(⌈/s←⍴y)×~l     the result shape
  t⊤t⍴⍳×/t              (1⌽⍳1+≢t)⍉↑⍳t, the index for each result element
  s⊥l                   the “del” in equation 2

(Misalignments in the display are due to defects in the APL Chat Forum software.)

For example:

      x←   2  1  2  0  1
      y← ? 5 13 19 17 11 ⍴ 100

      ⊢ l←x∘.=⍳0⌈1+⌈/x
0 0 1
0 1 0
0 0 1
1 0 0
0 1 0
      ⊢ t←s⌊.+(⌈/s←⍴y)×~l
17 11 5
      ⊢ s⊥l
11 3554 46376

      q← t⊤t⍴⍳×/t
      ⍴q
3 17 11 5
      q ≡ (1⌽⍳1+≢t)⍉↑⍳t
1

      (x⍉y) ≡ x tr79 y
1

Around this time (1981), Arthur Whitney was working with Ken Iverson on a model of APL written in APL [Iverson and Whitney 1982]. Thinking that perhaps he needed a model of transpose and secretly hoping that he might be impressed, I e-mailed him the 1-line transpose. His reply was very short, something like, it works.

1987

A model of ⍉ appeared in Some Uses of { and } [Hui 1987, §3.1]:

      ⍺⍉⍵:  ((>{⍳¨>(⍴⍵)⌊.+¯×~l)+.×(⍴⍵)⊥l){,⍵ ⊣ l←⍺∘.=⍳1+⌈/⍺

The function is a direct definition and is in “Dictionary APL” [Iverson 1987] wherein:

      {       a function
      {⍵      ⊃∘.,⌿⍵ (Cartesian product)
      ⍺{⍵     (⊂⍺)⌷⍵
      ¨       a dyadic operator (under AKA dual)
      ⍳¨>     ⍳¨ (iota each)
      >       ↑ (disclose AKA mix)
      ¯       ∞ (infinity)

tr16a in the 2016 section below is a close translation of this function in current Dyalog APL.

1992

From An Implementation of J [Hui 1992, p.60]:

      mask  =. =/ i.@>:@(>./)
      vec   =. >@{@:(i.&.>)@((<./ .+) _&*@-.)
      ind   =. vec +/ .* (#. |:)
      canta =. ($@] ind mask@[) { ,@]

2016

From A History of APL in 50 Functions [Hui 2016, §30].

Dyadic transpose, x⍉y, is probably one of the last primitives to be mastered for an APLer, but is actually straightforward to describe. Assign a distinct letter for each unique integer in x:

      0 1 2 3 …
      i j k l

If z←x⍉y, then z[i;j;k;…] equals y indexed by the letters corresponding to elements of x. For example:

      y← ? 5 13 19 17 11 ⍴ 100
      x←   2  1  2  0  1
      ⍝    k  j  k  i  j

      z←x⍉y
      i←?17 ⋄ j←?11 ⋄ k←?5
      z[i;j;k] = y[k;j;k;i;j]
1
      z[i;j;k] = y[⊂⍎¨'ijk'[x]]
1

From the above it can be seen that:

• the rank of z is 0⌈1+⌈/x
• the shape of z is (⍴y)⌊.+(⌈/⍴y)×x∘.≠⍳0⌈1+⌈/x [Berry 1979, p. 146]

And from that it is but a short step to a full model of x⍉y, essentially a problem of producing the ravelled representation from the strided representation [Nickolov 2013]. The version here is from [Hui 1987, §3.1], rendered in Dyalog APL:

      tr16a ← {(,⍵)⌷⍨⊂(↑,¨⍳(⍴⍵)⌊.+(⌈/⍴⍵)×~b)+.×(⍴⍵)⊥b←⍺∘.=⍳0⌈1+⌈/⍺}

The definition is briefer in an APL that supports infinity:

      tr16b ← {(,⍵)⌷⍨⊂(↑,¨⍳(⍴⍵)⌊.+∞×~b)+.×(⍴⍵)⊥b←⍺∘.=⍳0⌈1+⌈/⍺}

2017

From Tests, Derivations, Proofs [Hui 2017], presented at Functional Conf 2017, Bangalore, India, 2017-11-18. This version provides checks on the arguments and the result.

      ⍝
tr17 ← {⍵[(⊂⊂⍺) ⌷¨ ,¨ ⍳ ⍺[⍋⍺] {⌊/⍵}⌸ (⍴⍵)[⍋⍺]]}

assert←{⍺←'assertion failure' ⋄ 0∊⍵:⍺ ⎕signal 8 ⋄ shy←0}

TC←{
  assert 1=⍴⍴⍺ :
  assert (≢⍺)=⍴⍴⍵ :
  assert 0≤⍺ :                              ⍝ preconditions
  assert ⍺≡⌊⍺ :
  r←1+0⌈⌈/⍺
  assert ⍺(∊,∊⍨)⍳r :

  z←⍺ ⍺⍺ ⍵

  assert (⍴⍴z) = r :
  assert (⍴z)  ≡ ⍺[⍋⍺] {⌊/⍵}⌸ (⍴⍵)[⍋⍺] :    ⍝ post-conditions
  assert z[⊂v] = ⍵[⊂v[⍺]] ⊣ v←?⍴z :
  z
}

For example:

      x←   2  1  2  0  1
      y← ? 5 13 19 17 11 ⍴ 100

      (x⍉y) ≡ x ⍉TC y
1
      (x⍉y) ≡ x tr17 y
1
      (x⍉y) ≡ x tr17 TC y
1

2020

From APL Since 1978 [Hui and Kromberg 2020], §9.1 Array Representations.

APL arrays are stored in row-major (ravelled) order—consecutive items in a row are stored in consecutive memory locations; more generally, consecutive cells, and consecutive items in a cell, are in consecutive memory locations. Equivalently, in terms of major cells, array items are ordered in memory by major cells, by major cells within major cells, etc. Due to CPU caching on modern architectures, sequential memory access is more efficient than non-sequential access, and for good performance data design needs to take the order into account (favoring inverted tables [Hui 2018b], for example). In contrast, the Iliffe vector representation [Iliffe 1961] of a rank-n array (for n≥2) specifies a vector of pointers to rank-(n-1) subarrays (hello, major cells); as a result, data elements in different rows (cells) can be arbitrarily far from each other in memory. Iliffe vectors are implemented in Java, Perl, Python, and other programming languages.

Another representation is “strided representation” [Abrams 1970, §III; Munsey 1977; Nickolov 2013], where the ravelled data items are accessed through the function

      array[⊂⍵] ←→ (,array)[offset+⍵+.×stride]

Here, ⍵ are the usual indices i,j,k, … in each dimension, offset is an integer scalar, and stride is an integer vector with length the rank of array. This is a generalization of the usual index function where offset←0 and stride←×⍀⍢⊖1↓1,⍨⍴array or, equivalently,

      array[⊂⍵] ←→ (,array)[(⍴array)⊥⍵].

Strided representation permits terse descriptions of the selection functions take, drop, reverse, rotate, and transpose, which can be effected by computing the shape of the result and appropriate values for offset and stride. For example, the dreaded dyadic transpose ⍺⍉⍵ obtains as follows [Abrams 1970, p.73; Berry 1979, p.146; Hui 1987, §3.1; Hui 2016a, §30]:

      tr20←{
        shape  ← ⍺[⍋⍺] {⌊⌿⍵}⌸ (⍴⍵)[⍋⍺]
        stride ← ⍺[⍋⍺] {+⌿⍵}⌸ (×⍀⍢⊖1↓1,⍨⍴⍵)[⍋⍺]
        offset ← 0
        (,⍵)[offset+(↑,¨⍳shape)+.×stride]
      }

The last line of tr20 produces the row-major representation from the strided representation. The computation of the strided representation itself is O(≢⍺), or O(≢⍴⍵) since (≢⍺)=≢⍴⍵. …

(The expression ×⍀⍢⊖1↓1,⍨⍴⍵ makes use of the as yet unimplemented under operator, and is equivalent to ⊖×⍀⊖1↓1,⍨⍴⍵.)

A Transpose Puzzle

On 2003-08-03 Kip Murray posed a puzzle on the J Programming Forum. The following is a slight rewording of it:

An array A is symmetric if A≡p⍉A for any permutation p of order r←≢⍴A. To check whether an array is symmetric, can you do better than checking every permutation?

“Symmetric Array” is a generalization of symmetric matrix: A matrix M is symmetric if M≡⍉M. But ⍉M ↔ 1 0⍉M and of course M≡0 1⍉M, so M is symmetric if M≡p⍉M for any permutation p of order 2.

An example of a symmetric array:

      A← ⊃ ∘.+⌿ 5⍴⊂?5⍴100
      ⍴A
5 5 5 5 5
      A ≡ (5?5)⍉A
1

You can check for symmetry by checking every permutation. (Function perm is from the Dyalog Blog post Permutations, 2015-07-20.)

      perm←{0=⍵:1 0⍴0 ⋄ ,[⍳2](⊂0,1+∇ ¯1+⍵)⌷⍤1 ⍒⍤1 ∘.=⍨ ⍳⍵}

      ⊢ r← ≢⍴A
5
      ⍴ perm r
120 5
      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 ⊢A
120

Spoiler alert.
.
.
.
.
.
.
.
.
.
.

0. If p and q are permutations of order r, then

      (p⍉q⍉A) ≡ p[q]⍉A

1. Let P be a set of permutations of order r. From 0, the following are equivalent:

      P            (⊢ ≡ ⍉)⍤1 99 ⊢ A
      (subgroup P) (⊢ ≡ ⍉)⍤1 99 ⊢ A

where

      subgroup←{{∪(2 1*⍨⍴⍵)⍴⍵[;⍵]}⍣≡(⍉⍪⍳⊃⊖⍴⍵)⍪⍵}

produces the subgroup generated by permutation(s) ⍵. (Compare this with an alternative formulation, the operator C in the Limit and Closure post.)

2. The two permutations p←1 0,2↓⍳r (swapping the first two items, for r≥2) and q←1⊖⍳r (rotating by 1) generate the entire permutation group: p and q suffice to generate p1←p and q1←(1⌽⍳r-1),r-1, which are like p and q except ¯1↑⍵ is invariant under ⍵[p1] and ⍵[q1] [Hui 1987, §4.4].

The following examples illustrate the point:

      gen←{(2,⍵)⍴(1⊖⍳⍵),(1<⍵)⌿1 0,2↓⍳⍵}
      gen 5
1 2 3 4 0
1 0 2 3 4
      ⍴subgroup gen 5
120 5
      ⍴∘subgroup∘gen¨ ⍳8
┌───┬───┬───┬───┬────┬─────┬─────┬──────┐
│1 0│1 1│2 2│6 3│24 4│120 5│720 6│5040 7│
└───┴───┴───┴───┴────┴─────┴─────┴──────┘
      !⍳8
1 1 2 6 24 120 720 5040

3. Therefore, A is symmetric if both of the following are true:

      A ≡ A ⍉⍨ 1⊖⍳r←≢⍴A
      A ≡ A ⍉⍨ 1 0,2↓⍳r  ⍝ if r≥2

That is, it suffices to check two permutations instead of !≢⍴A permutations.

For example:

      ⊢ x←?5⍴1000
925 845 368 482 937
      A← x∘.+x∘.+x∘.+x∘.+x

      A ≡ A ⍉⍨ 1⊖⍳r←≢⍴A
1
      A ≡ A ⍉⍨ 1 0,2↓⍳r
1
      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 ⊢A
120

      ○@0⊢x
2905.97 845 368 482 937
      B← x∘.+x∘.+x∘.+(○@0⊢x)∘.+x

      B ≡ B ⍉⍨ 1⊖⍳r←≢⍴B
0
      B ≡ B ⍉⍨ 1 0,2↓⍳r
1
      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 ⊢B
24

      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 -@(⊂0 1 2 3 4)⊢A
1
      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 -@(⊂0 0 2 3 4)⊢A
2
      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 -@(⊂0 0 0 3 4)⊢A
6
      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 -@(⊂0 0 0 0 4)⊢A
24
      +⌿ (perm r) (⊢ ≡ ⍉)⍤1 99 -@(⊂0 0 0 0 0)⊢A
120

[Bob Armstrong]

Roger , I greatly enjoy your histories of various APL concepts .
I remember being quite proud when I got a 4 dimensional transpose to work . It always seemed like the index sequence needed was sort of inverse to the way I thought . But I could see the reasons . I think I may have also managed to get some particular trace to work .

I know I expressed long ago that the greatest question I had when I followed Arthur Whitney's path thru K was : does 1st axis and simple ' flip ( K : monadic + ) between adjacent ` dimensions ( depths ) suffice ?

CoSy is much more concerned with the simplest possible implementation of structures ( list of lists ala Whitney ( are they equivalent to Iliffe ? ) where the only structures are homogeneous lists , and , essentially , lists of pointers ) and essential vocabulary . Thus no deeper calculable memory structure .

Given Forth's simple white space is the prime delimiter and any non-blank string being a valid name , I'll wait for a unicode interface so people can use any symbol sets for names they like .

But , there's nothing to prevent one from implementing fancier concepts .

Bottom line , I'd say , yes the simple structure and ' flip is sufficient .

A somewhat different issue comes up given CoSy's modulo indexing : What's the flip of a ragged list ?

Here's the current definition of ' flip :

Code: Select all

: flip ( CSob -- CSob )    | Transpose list of 2 lists .
| returns list of each item of 1th list w corresponding item of 1st suject
| to the minimum length of the 2 lists .
   dup @ if ;then      | transpose of a simple obj is itself
   dup i# 0;drop       | same for empty
   refs+> dup ' rho 'm ,/ ' min _./ >_ cellVecInit >aux>    | ob nbdy
    i# 0 ?do dup i _nth refs+> aux@ i i! loop
    refs- aux> ;

where

Code: Select all

: _nth ( CSob n -- CSadr )     | extracts nth item from each item of CSob 
   _i 2refs+> ev  2 pick i# 0 ?do 2 pick i i@ 2 pick at\ cL loop >r 2refs- r> ;
 

My main question is whether to use the minimum or maximum of the counts of the lists being flipped . Most operations in CoSy use maximum , using modulo indexing to naturally extent the shorter to the longer . But it took me trying that to see it's generally useless .

But the need to flip ragged lists rather than paired lists is at most rare .