A user recently asked how to find the set of indices corresponding to a given value of a two-dimensional Array.

There are several ways to handle this problem.  For  a single value, a simple scan through the Array suffices:

FindIndices1 := proc(A :: Array, val)
local i,j,irng,jrng;
    (irng,jrng) := rtable_dims(A);
    {seq(seq(`if`(A[i,j]=val, [i,j], NULL), i=irng), j=jrng)};
end proc:

With a multi-core machine, the ?Threads package may be used to perform the search in parallel:

FindIndices2 := proc(A :: Array, val)
local i,j,irng,jrng;
    (irng,jrng) := rtable_dims(A);
    {Threads:-Seq(SearchRow(A,i,jrng,val),i=irng)};
end proc:

SearchRow := proc(A, i, jrng, val)
local j;
    seq(`if`(A[i,j] = val, [i,j], NULL), j=jrng);
end proc:

Another approach is to use ?rtable_scanblocks. That has the advantage of working with Arrays of any dimension:

FindIndices3 := proc(A :: rtable, val)
local T;
    T := table();
    rtable_scanblock(A, [], 'passindex'
                     , proc(x,indx)
                           if x=val then
                               T[indx] := NULL;
                           end if;
                       end proc
                    );
    {indices(T,'nolist')};
end proc:

The procedures just described are suitable for searching for a few values. If you have to look up many values, they are not efficient in that each search requires a complete scan of the Array.  A better approach is to generate a table that maps each value in the Array to the set of indices to which it corresponds.   One way, reasonable if the Array has few duplicate entries, is to scan the Array and append each element to a sequence in a table.

InverseLookup1 := proc(A::Array)
local irng, jrng, i, j, cnt, val, valToIndex, cntToVal;
    # Check Array
    if not rtable_num_dims(A) = 2 then
        error "Array must be two-dimensional";
    end if;  
    (irng,jrng) := rtable_dims(A);
    if lhs(irng)<>1 or lhs(jrng)<>1 then
        error "Array indices must start at 1";
    end if;
    # Append indices to a sequence in a table
    cnt := 0;
    for i to rhs(irng) do
        for j to rhs(jrng) do
            val := A[i,j];
            if assigned(valToIndex[val]) then
                valToIndex[val] := valToIndex[val], [i,j];
            else
                cnt := cnt+1;
                cntToVal[cnt] := val;
                valToIndex[val] := [i,j];
            end if;
        end do;
    end do;
    # Convert sequences to sets
    for cnt to cnt do
        val := cntToVal[cnt];
        valToIndex[val] := { valToIndex[val] };
    end do;
    return eval(valToIndex);
end proc:

A problem with this technique is that its worst-case performance, for an m x n Array, is O(m^2*n^2), which occurs when all the values are identical.  This can be seen with

CodeTools:-Usage(InverseLookup1(Array(1..1000,1..4,'fill'=1))):
memory used=63.26MiB, alloc change=40.74MiB, cpu time=0.32s, real time=0.32s
CodeTools:-Usage(InverseLookup1(Array(1..2000,1..4,'fill'=1))):
memory used=246.01MiB, alloc change=58.86MiB, cpu time=1.12s, real time=1.12s

This behavior comes from the incremental assignment to valToIndex in the nested do-loop.  The usual way to break this square-law performance is to assign distinct quantities (here indices), to separate entries in a table, and then combine them after all have been entered.  There is, however, a related but more efficient technique that can be used.

Consider a two-dimensional Array with the value 1 stored at indices [1,1], [2,2], [3,3], [4,4].   Assume this is the order in which they are scanned (other indices are also scanned, but they hold different values).   When 1 is found at [1,1], we create a link

1 --> [1,1]

When the second 1 is scanned at [2,2] we insert this at the start of the chain

1 --> [2,2] --> [1,1]

When finished scanning we have the chain

1 --> [4,4] --> [3,3] --> [2,2] --> [1,1]

To implement linked lists in Maple, we use a table, valToHead, that maps a value to the first index.  The remaining indices in the linked list are stored in a three-dimensional Array, where the entries for ([i,j,1], [i,j,2]) are the indices for the next link.  The end of the chain is denoted by a 0 for the [i,j,1] entry.

InverseLookup2 := proc(A::Array)
local irng,jrng,i,j,k,val,cnt,indx, Links, valToHead, valToIndx, CntToVal, T;
    # Check Array
    if not rtable_num_dims(A) = 2 then
        error "Array must be two-dimensional";
    end if;
    (irng,jrng) := rtable_dims(A);
    if lhs(irng)<>1 or lhs(jrng)<>1 then
     error "Array indices must start at 1";
    end if;
    # Initialize Links Array
    Links  := Array(irng,jrng,1..2, 'datatype'=integer[4]);
    valToHead := table();
    cnt := 0;
    # Create links
    for i to rhs(irng) do
        for j to rhs(jrng) do
            val := A[i,j];
            if assigned(valToHead[val]) then
                (Links[i,j,1],Links[i,j,2]) := valToHead[val];
            else
                cnt := cnt+1;
                CntToVal[cnt] := val;
            end if;
            valToHead[val] := (i,j);
        end do;
    end do;
    # Convert links to a table of sets
    for cnt to cnt do
        val := CntToVal[cnt];
        indx := valToHead[val];
        for k do
            T[k] := [indx];
            indx := (Links[indx,1],Links[indx,2]);
            if indx[1] = 0 then break; end if;
        end do;
        valToIndx[val] := { seq(T[k], k=1..k) };
    end do;
    return eval(valToIndx);
end proc:

The performance advantage of this technique compared to the worst-case performance of InverseLookup1 is substantial:

CodeTools:-Usage(InverseLookup2(Array(1..1000,1..4,'fill'=1))):
memory used=2.06MiB, alloc change=0 bytes, cpu time=0.05s, real time=0.04s
CodeTools:-Usage(InverseLookup2(Array(1..2000,1..4,'fill'=1))):
memory used=4.20MiB, alloc change=0 bytes, cpu time=0.05s, real time=0.05s

At the cost of some speed, we can use rtable_scanblock to implement this same technique so that it works with an Array of any dimension.

InverseLookup3 := proc(A::Array)
local dims,k,val,cnt,indx,ndims, Links, valToHead, valToIndx, cntToVal, T;
    dims := [rtable_dims(A)];
    if not andmap(rng -> lhs(rng)=1, dims) then
        error "Array indices must start at 1";
    end if;
    ndims := nops(dims);
    Links  := Array(dims[], 1..ndims, 'datatype'=integer[4]);
    valToHead := table();
    cnt := 0;
    # Create the links
    rtable_scanblock(A, [], 'passindex'
                     , proc(val,indx)
                       local i,head;
                           if assigned(valToHead[val]) then
                               head := valToHead[val];
                               for i to ndims do
                                   Links[indx[],i] := head[i];
                               end do;
                           else
                               cnt := cnt+1;
                               cntToVal[cnt] := val;
                           end if;
                           valToHead[val] := indx;
                       end proc
                    );
    # Create the output table by following the links for each value.
    for cnt to cnt do
        val := cntToVal[cnt];
        indx := valToHead[val];
        for k do
            T[k] := indx;
            indx := [seq(Links[indx[],i], i=1..ndims)];
            if indx[1] = 0 then break; end if;
        end do;
        valToIndx[val] := { seq(T[k], k=1..k) };
    end do;
    return eval(valToIndx);
end proc:

Please Wait...