(* Copyright 1999-2001, Patrick Hunziker

  This library is free software; you can redistribute it and/or modify it under the terms of the
  GNU Library General Public License as published by the Free Software Foundation;
  either version 2 of the License, or any later version.

  This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  See the GNU Library General Public License for more details.

  You should have received a copy of the GNU Library General Public License along with this library;
  if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

Patrick Hunziker,Basel.
email Patrick.Hunziker@unibas.ch
*)

(** Version 1.0, 19.1.2001 *)

MODULE MultiArrayRiders;  (** Patrick Hunziker, Basel,  **)
(** Implements an array rider access mechanism for multidimensional arrays of arbitrary
    dimensions defined in MultiArrays*)
IMPORT MultiArrays, Out:= Console, Input := Kernel;
CONST (** behaviour of array rider at end of array line;
        not yet completely implemented.
        The seemingly more exotic variants are especially useful in image processing *)
  halt = 0;
  zeropadding = 1;  (* not yet implemented *)
  constant = 2;    (* not yet implemented *)
  circular* = 3;  (** after finishing one line, the same line is restarted *)
  mirror = 4;      (* not yet implemented *)
  incremental* = 5; (** after finishing one line, the next line is started *)

  (** rider has not passed any border of the array *)
  noteol* = MAX(LONGINT);

TYPE
    (** Array riders allow traversal of arbitrary dimensional array using procedures Inc() and Dec(),
        and can be positioned using Set(). *)
    Rider* = RECORD
      array-: MultiArrays.Array; (** points to the array the rider is based on *)
      order-: LONGINT;   (** dimensionality of array *)
      offset- : LONGINT;  (** Rider position in linear array representation *)
      eol*: LONGINT;     (** Rider has gone beyond the border of the
                       line of indicated dimension .
                       if eol=noteol, rider is inside array *)
      eolBehaviour*: INTEGER; (** What to do when reaching eol. Not yet completely satisfactory *)
      dimension,        (* dimensions of Array, in vector notation *)
      pos,              (* position of rider in Array, in vector notation *)
      step:     MultiArrays.SizeVector;  (* unit increment for offset in each dimension, in vector notation *)
      END;

    PROCEDURE CalculatePos (pos, dimension: MultiArrays.SizeVector): LONGINT;
      VAR maxI, res, i: LONGINT;
      BEGIN
        maxI := LEN(dimension^)-1;
        ASSERT(LEN(pos^) = LEN(dimension^));
        res := pos[maxI];
        FOR i := 1 TO maxI DO res := res*dimension[maxI-i]+pos[maxI-i] END;
        RETURN res
      END CalculatePos;

    PROCEDURE InitRider* (VAR R: Rider; A: MultiArrays.Array; pos: MultiArrays.SizeVector);
            (** Sets array rider in position pos in array A *)
      VAR i, step: LONGINT;
      BEGIN
        ASSERT(MultiArrays.Order(A) = LEN(pos^));
        R.array := A;
        R.order := MultiArrays.Order(A);
        NEW(R.pos,R.order);
        NEW(R.step,R.order);
        NEW(R.dimension,R.order);
        step := 1;
        FOR i := 0 TO R.order-1 DO
          ASSERT(pos[i] <= MultiArrays.Len(A,i));
          R.pos[i] := pos[i];
          R.step[i] := step; step := step*MultiArrays.Len(A,i);
          R.dimension[i] := MultiArrays.Len(A,i)
        END;
        R.eol := noteol;
        R.offset := CalculatePos(R.pos,MultiArrays.Size(A));
        R.eolBehaviour := incremental
      END InitRider;

    PROCEDURE SetRider* (VAR R: Rider; pos: MultiArrays.SizeVector);
    VAR i: LONGINT;
    BEGIN
      ASSERT(R.array # NIL);
      ASSERT(LEN(pos^) = R.order);
      FOR i := 0 TO R.order-1 DO ASSERT(pos[i] < R.dimension[i]); R.pos[i] := pos[i] END;
      R.offset := CalculatePos(pos,R.dimension);
      R.eol := noteol
    END SetRider;

    PROCEDURE Inc* (VAR R: Rider; Dim: LONGINT);
            (** array rider advances one element in dimension Dim;
        at end of line, eol is assigned the number of the dimension overflown *)
      BEGIN
        ASSERT(Dim < R.order);
        IF R.pos[Dim] < R.dimension[Dim]-1
        THEN INC(R.pos[Dim]); INC(R.offset, R.step[Dim]);
        ELSE
          R.eol := Dim;
          CASE R.eolBehaviour OF
            halt: HALT(100);
            | zeropadding: HALT(100); (* not yet implemented *)
            | constant:
            | mirror: HALT(100); (* not yet implemented *)
            | incremental:
              R.pos[Dim] := 0;
              IF Dim < R.order-1
              THEN
                INC(R.offset, R.step[Dim]-R.step[Dim+1]);
                Inc(R, R.eol+1)
              ELSE INC(R.offset, R.step[Dim]-R.array.len)
              END
            | circular:
              R.pos[Dim] := 0;
              IF Dim < R.order-1
              THEN INC(R.offset, R.step[Dim]-R.step[Dim+1])
              ELSE INC(R.offset, R.step[Dim]-R.array.len)
              END
            ELSE HALT(100)
          END
        END
      END Inc;

    PROCEDURE Dec* (VAR R: Rider; Dim: LONGINT);
            (** array rider goes back one element in dimension Dim *)
      BEGIN
        ASSERT(Dim < R.order);
        IF R.pos[Dim] > 0
        THEN DEC(R.pos[Dim]); DEC(R.offset, R.step[Dim]);
        ELSE R.eol := Dim;
          CASE R.eolBehaviour OF
            halt: HALT(100);
            | zeropadding: HALT(100); (* not yet implemented *)
            | constant:
            | mirror: HALT(100); (* not yet implemented *)
            | incremental:
              R.pos[Dim] := R.dimension[Dim]-1;
              IF Dim > 0
              THEN
                DEC(R.offset, R.step[Dim]-R.step[Dim+1]);
                Dec(R, R.eol+1)
              ELSE DEC(R.offset, R.step[Dim]-R.array.len)
              END
            | circular:
              R.pos[Dim] := R.dimension[Dim]-1;
              IF Dim > 0
              THEN DEC(R.offset, R.step[Dim]-R.step[Dim+1])
              ELSE DEC(R.offset, R.step[Dim]-R.array.len)
              END
            ELSE HALT(100)
          END
        END
      END Dec;

    PROCEDURE Pos* (R: Rider): MultiArrays.SizeVector;
            (** gives actual position of R in its associated array *)
      VAR i: LONGINT; res: MultiArrays.SizeVector;
      BEGIN
        NEW(res,R.order);
        FOR i := 0 TO R.order-1 DO res[i] := R.pos[i] END;
        RETURN res
      END Pos;

(** elementwise reading from Array Rider, followed by advancing the rider by
one step in direction "dir"; with specific "eolBehaviour" (see above) at border of array *)
    PROCEDURE ReadSInt* (VAR R: Rider; dir: LONGINT; VAR s: SHORTINT);
    BEGIN
      IF R.array IS MultiArrays.SIntArray
      THEN s := R.array(MultiArrays.SIntArray).s[R.offset]; Inc(R, dir)
      ELSE HALT(100)
      END
    END ReadSInt;

    PROCEDURE ReadSIntRun* (VAR R: Rider; dir: LONGINT;
                              VAR srun: ARRAY OF SHORTINT; n: LONGINT);
    VAR i, step, offset, pos, dim: LONGINT; array: MultiArrays.SIntArray;
    BEGIN
      ASSERT(LEN(srun) >= n);
      ASSERT(dir < R.order);
      ASSERT(R.array IS MultiArrays.SIntArray);
      array := R.array(MultiArrays.SIntArray);
      offset := R.offset;
      step := R.step[dir];
      pos := R.pos[dir];
      dim := R.dimension[dir];
      CASE R.eolBehaviour OF
         halt: HALT(100); (* not yet implemented *)
      | incremental:
        IF offset+(n-1)*step < R.array.len
        THEN FOR i := 0 TO n-1 DO srun[i] := array.s[offset]; INC(offset, step) END
        ELSE HALT(100) (* not yet implemented *)
        END
      | circular:
        IF R.pos[dir]+n-1 < dim
        THEN FOR i := 0 TO n-1 DO srun[i] := array.s[offset]; INC(offset, step) END
        ELSE
          FOR i := 0 TO n-1 DO
            srun[i] := array.s[offset+((pos+i) MOD dim)*step] (* can further be optimized *)
          END
        END
      ELSE HALT(100) (* not yet implemented *)
      END
    END ReadSIntRun;

    PROCEDURE ReadInt* (VAR R: Rider; dir: LONGINT; VAR i: INTEGER);
    BEGIN
      IF R.array IS MultiArrays.IntArray
      THEN i := R.array(MultiArrays.IntArray).i[R.offset]; Inc(R, dir);
      ELSE HALT(100)
      END
    END ReadInt;

    PROCEDURE ReadLInt* (VAR R: Rider; dir: LONGINT; VAR j: LONGINT);
    BEGIN
      IF R.array IS MultiArrays.LIntArray
      THEN j := R.array(MultiArrays.LIntArray).j[R.offset]; Inc(R, dir)
      ELSE HALT(100)
      END
    END ReadLInt;

(*    PROCEDURE ReadHInt* (VAR R: Rider; dir: LONGINT; VAR h: HUGEINT);
    BEGIN
      HALT(100) (* yet to implement *)
    END ReadHInt; *)

    PROCEDURE ReadReal* (VAR R: Rider; dir: LONGINT; VAR x: REAL);
    BEGIN
      IF R.array IS MultiArrays.RealArray
      THEN x := R.array(MultiArrays.RealArray).x[R.offset]; Inc(R, dir)
      ELSE HALT(100)
      END;
    END ReadReal;

    PROCEDURE ReadRealRun* (VAR R: Rider; dir: LONGINT;
                              VAR srun: ARRAY OF REAL; n: LONGINT);
    VAR i, step, offset, pos, dim: LONGINT; array: MultiArrays.RealArray;
    BEGIN
      ASSERT(LEN(srun) >= n);
      ASSERT(dir < R.order);
      ASSERT(R.array IS MultiArrays.RealArray);
      array := R.array(MultiArrays.RealArray);
      offset := R.offset;
      step := R.step[dir];
      pos := R.pos[dir];
      dim := R.dimension[dir];
      CASE R.eolBehaviour OF
         halt: HALT(100); (* not yet implemented *)
      | incremental:
        IF offset+(n-1)*step < R.array.len
        THEN FOR i := 0 TO n-1 DO srun[i] := array.x[offset]; INC(offset, step) END
        ELSE HALT(100) (* not yet implemented *)
        END;
      | circular:
        IF R.pos[dir]+n-1 < dim
        THEN FOR i := 0 TO n-1 DO srun[i] := array.x[offset]; INC(offset, step) END
        ELSE
          FOR i := 0 TO n-1 DO
            srun[i] := array.x[offset+((pos+i) MOD dim)*step] (* can further be optimized *)
          END
        END
      ELSE HALT(100) (* not yet implemented *)
      END
    END ReadRealRun;

    PROCEDURE ReadLReal* (VAR R: Rider; dir: LONGINT; VAR y: LONGREAL);
    BEGIN
      IF R.array IS MultiArrays.LRealArray
      THEN y := R.array(MultiArrays.LRealArray).y[R.offset]; Inc(R, dir)
      ELSE HALT(100)
      END
    END ReadLReal;

(*    PROCEDURE ReadBool* (VAR R: Rider; dir: LONGINT; VAR b: BOOLEAN);
    BEGIN
      HALT(100) (* to implement *)
    END ReadBool; *)

(*    PROCEDURE ReadComplex* (VAR R: Rider; dir: LONGINT; VAR z: MultiArrays.Complex);
    BEGIN
      HALT(100) (* yet to implement *)
    END ReadComplex; *)

    PROCEDURE WriteSInt* (VAR R: Rider; dir: LONGINT; s: SHORTINT);
    BEGIN
      IF R.array IS MultiArrays.SIntArray
      THEN R.array(MultiArrays.SIntArray).s[R.offset] := s; Inc(R, dir)
      ELSE HALT(100) END
    END WriteSInt;

    PROCEDURE WriteInt* (VAR R: Rider; dir: LONGINT; i: INTEGER);
    BEGIN
      IF R.array IS MultiArrays.IntArray
      THEN R.array(MultiArrays.IntArray).i[R.offset] := i; Inc(R, dir)
      ELSE HALT(100) END
    END WriteInt;

    PROCEDURE WriteLInt* (VAR R: Rider; dir: LONGINT; j: LONGINT);
    BEGIN
      IF R.array IS MultiArrays.LIntArray
      THEN R.array(MultiArrays.LIntArray).j[R.offset] := j; Inc(R, dir)
      ELSE HALT(100) END
    END WriteLInt;

(*    PROCEDURE WriteHInt* (VAR R: Rider; dir: LONGINT; h: HUGEINT);
    BEGIN
      HALT(100); (* yet to implement *) END
    END WriteHInt; *)

    PROCEDURE WriteReal* (VAR R: Rider; dir: LONGINT; x: REAL);
    BEGIN
      IF R.array IS MultiArrays.RealArray
      THEN R.array(MultiArrays.RealArray).x[R.offset] := x; Inc(R, dir)
      ELSE HALT(100) END
    END WriteReal;

    PROCEDURE WriteRealRun* (VAR R: Rider; dir: LONGINT; srun: ARRAY OF REAL; n: LONGINT);
    VAR i, step, offset, pos, dim: LONGINT; array: MultiArrays.RealArray;
    BEGIN
      ASSERT(LEN(srun) >= n);
      ASSERT(dir < R.order);
      ASSERT(R.array IS MultiArrays.RealArray);
      array := R.array(MultiArrays.RealArray);
      offset := R.offset;
      step := R.step[dir];
      pos := R.pos[dir];
      dim := R.dimension[dir];
      CASE R.eolBehaviour OF
         halt: HALT(100); (* not yet implemented *)
      | incremental:
        IF offset+(n-1)*step < R.array.len
        THEN
          FOR i := 0 TO n-1 DO array.x[offset] := srun[i]; INC(offset, step) END
        ELSE HALT(100) (* not yet implemented *)
        END
      | circular:
        IF R.pos[dir]+n-1 < dim
        THEN
          FOR i := 0 TO n-1 DO array.x[offset] := srun[i]; INC(offset, step) END
        ELSE
          FOR i := 0 TO n-1 DO
            array.x[offset+((pos+i) MOD dim)*step] := srun[i] (* can further be optimized *)
          END
        END
      ELSE HALT(100) (* not yet implemented *)
      END
    END WriteRealRun;

    PROCEDURE WriteLReal* (VAR R: Rider; dir: LONGINT; y: LONGREAL);
    BEGIN
      IF R.array IS MultiArrays.LRealArray
      THEN R.array(MultiArrays.LRealArray).y[R.offset] := y; Inc(R, dir)
      ELSE HALT(100) END
    END WriteLReal;

    PROCEDURE WriteBool* (VAR R: Rider; dir: LONGINT; b: BOOLEAN);
    BEGIN
      IF R.array IS MultiArrays.BoolArray
      THEN R.array(MultiArrays.BoolArray).b[R.offset] := b; Inc(R, dir)
      ELSE HALT(100) END
    END WriteBool;

(*    PROCEDURE WriteComplex* (VAR R: Rider; dir: LONGINT; VAR z: MultiArrays.Complex);
    BEGIN
      HALT(100) (* yet to implement *)
    END WriteComplex; *)

    PROCEDURE InvertSign (s: SHORTINT): SHORTINT; (* Test procedure for unary operations *)
    BEGIN
      RETURN -s
    END InvertSign;

    PROCEDURE Assign (VAR S: SHORTINT;s: SHORTINT); (* Testing *)
    BEGIN
      S := s
    END Assign;


    PROCEDURE Test*; (** Tests if eol mechanism is working correctly *)
      VAR pos, dimension: MultiArrays.SizeVector; A: MultiArrays.Array;
      i, j: LONGINT; R: Rider;
      BEGIN
        MultiArrays.SizeVector4(dimension, 10, 10, 10, 10);
        MultiArrays.SizeVector4(pos, 2, 3, 4, 5);
        NEW(A); MultiArrays.InitInt(A, dimension, NIL, FALSE);
        InitRider(R,A,pos);
        R.eolBehaviour := circular;
        FOR j := 0 TO 3 DO
          FOR i := 1 TO 10 DO
            Inc(R, j); Out.Int(CalculatePos(R.pos,MultiArrays.Size(A)), 5); Out.Ln;
            IF R.eol#noteol THEN
              Out.String("R.eol:"); Out.Int(R.eol, 5); Out.Ln;
              R.eol := noteol
            END
          END;
          Out.String("----"); Out.Ln
        END
      END Test;

      PROCEDURE Test1*;
      VAR A1: MultiArrays.Array;
          SA: ARRAY 256 OF SHORTINT;
          R1: Rider;
          dim0, dim1: MultiArrays.SizeVector;
          starttime, endtime, a, b, c, d, opcount1, opcount2: LONGINT;
          A3: POINTER TO ARRAY OF ARRAY OF ARRAY OF ARRAY OF SHORTINT;
          s: SHORTINT;
      BEGIN
        Out.Ln;
        Out.String("**********************************"); Out.Ln;
        Out.String(" Benchmark:"); Out.Ln;
        Out.String(" Arbitrary arrays with riders vs. ARRAY[x,y,z,...] concept"); Out.Ln;
        Out.String("----------------------------------"); Out.Ln;
        Out.String("----------------------------------"); Out.Ln;
        MultiArrays.SizeVector4(dim0, 0, 0, 0, 0);
        MultiArrays.SizeVector4(dim1, 256, 128, 8, 8);
        MultiArrays.InitSInt(A1, dim1, NIL, FALSE);
        InitRider(R1, A1, dim0);
        R1.eolBehaviour := incremental;
        opcount1 := 0;
(* ASSIGN *)
        starttime := Input.Time();
        REPEAT          (* main loop of elementwise rider writing *)
              ReadSInt(R1,0,s);
              INC(opcount1);
            UNTIL R1.eol=R1.order-1;
        endtime := Input.Time();
        Out.String("ASSIGN:");Out.Ln;
        Out.String("Arbitrary multidimensional array: elementwise writing:"); Out.Ln;
        Out.String("  time: "); Out.Int(endtime-starttime,5);
        Out.String("  opcount: "); Out.Int(opcount1, 5); Out.Ln;
        Out.String("----------------------------------"); Out.Ln;
        MultiArrays.SizeVector4(dim0,0,0,0,0);
        MultiArrays.SizeVector4(dim1,256,128,8,8);
        SetRider(R1,dim0); R1.eolBehaviour := circular;
        opcount2 := 0;
        dim0[0] := 0;
        starttime := Input.Time();
            FOR b := 0 TO dim1[1]-1 DO  (* main loop of linewise reading *)
              dim0[1] := b;
              FOR c := 0 TO dim1[2]-1 DO
                dim0[2] := c;
                FOR d := 0 TO dim1[3]-1 DO
                  dim0[3] := d;
                  SetRider(R1,dim0);
                  ReadSIntRun(R1,0,SA,256);
                  INC(opcount2, 256)
                END
              END
            END;
        endtime := Input.Time();
        Out.String(" Arbitrary multidimensional array: line reading with rider"); Out.Ln;
        Out.String("  time:  ");
        Out.Int(endtime-starttime, 5);
        Out.String("  opcount: "); Out.Int(opcount2, 5); Out.Ln;
        Out.String("----------------------------------"); Out.Ln;

        NEW(A3, 256, 128, 8, 8 );
        starttime := Input.Time();
        opcount2 := 0;
        FOR a := 0 TO LEN(A3^,0)-1 DO   (* main loop of conventional array handling *)
              FOR b := 0 TO LEN(A3^,1)-1 DO
                FOR c := 0 TO LEN(A3^,2)-1 DO
                  FOR d := 0 TO LEN(A3^,3)-1 DO
                      SA[a] := A3[a,b,c,d]; INC(opcount2)
                    END
                END
              END
            END;
        endtime := Input.Time();
        Out.String(" conventional multidimensional array: index line reading:"); Out.Ln;
        Out.String("  time:  ");
        Out.Int(endtime-starttime, 5);
        Out.String("  opcount: "); Out.Int(opcount2, 5); Out.Ln;
        Out.String("**********************************"); Out.Ln
  END Test1;

(* Intel may have register problems with the following procedure *)
(*    PROCEDURE Test2*;
      VAR A1, A2: MultiArrays.Array;
          R1: Rider;
          dim0, dim1: MultiArrays.SizeVector;
          i, starttime, endtime, a, b, c, d, e, f, dir, opcount1, opcount2: LONGINT;
          A3: POINTER TO ARRAY OF ARRAY OF ARRAY OF ARRAY OF ARRAY OF ARRAY OF SHORTINT;
          s: SHORTINT;
      BEGIN
        Out.Ln;
        Out.String("**********************************"); Out.Ln;
        Out.String(" Benchmark:"); Out.Ln;
        Out.String(" Arbitrary arrays with riders vs. ARRAY[x,y,z,...] concept"); Out.Ln;
        Out.String("----------------------------------"); Out.Ln;
        NEW(dim0, 6); FOR i := 0 TO 5 DO dim0[i] := 0 END;
        NEW(dim1, 6);
        dim1[0] :=  64; dim1[1] :=  32; dim1[2] :=  16; dim1[3] :=  16;
        dim1[4] :=  2; dim1[5] :=  2;
        MultiArrays.InitSInt(A1, dim1, NIL, FALSE);
        InitRider(R1, A1, dim0);
        R1.eolBehaviour := incremental;
        opcount1 := 0;
(* ASSIGN *)
        starttime := Input.Time();
        REPEAT
          WriteSInt(R1,0,1);
          INC(opcount1);
        UNTIL R1.eol=5;
        endtime := Input.Time();
        Out.String("ASSIGN:");Out.Ln;
        Out.String("arbitrary array rider assignement:"); Out.Ln;
        Out.String("  time: "); Out.Int(endtime-starttime, 5);
        Out.String("  opcount: "); Out.Int(opcount1, 5); Out.Ln;
        Out.String("----------------------------------"); Out.Ln;

        NEW(A3, 64,32,16,16, 2,2);
        starttime := Input.Time();
        opcount2 := 0;
        FOR a := 0 TO LEN(A3^,0)-1 DO
          FOR b := 0 TO LEN(A3^,1)-1 DO
            FOR c := 0 TO LEN(A3^,2)-1 DO
              FOR d := 0 TO LEN(A3^,3)-1 DO
                FOR e := 0 TO LEN(A3^,4)-1 DO
                  FOR f := 0 TO LEN(A3^,5)-1 DO
                    Assign(A3[a,b,c,d,e,f],1); INC(opcount2)
                  END
                END
              END
            END
          END
        END;
        endtime := Input.Time();
        Out.String("multidim index assignement:"); Out.Ln;
        Out.String("  time:  ");
        Out.Int(endtime-starttime, 5);
        Out.String("  opcount: "); Out.Int(opcount2, 5); Out.Ln
END Test2;
*)

BEGIN
END MultiArrayRiders.

MultiArrayRiders.Test1 (4D), MultiArrayRiders.Test2 (6D)

Compares execution times for
                    -arbitrary dimensional approach with riders
                    -conventional ARRAY [x,y,z,...] approach
This is done for elementwise assignements and reading of "runs" of data.

