Re: [boost] [multi_array] Is there interest for an [index_list] library?

Subject: Re: [boost] [multi_array] Is there interest for an [index_list] library?
From: Larry Evans (cppljevans_at_[hidden])
Date: 2011-04-22 16:01:09

On 04/11/11 15:15, Pierre-Andre Noel wrote:
[snip]
>
> The "Domain" concept defined in the IndexList library generalizes the
> conversion you are talking about. Passing a number in the range
> [0,num_elements()-1] to the "expand()" member of a Domain returns a
> collection of indices, and passing that collection of indices to its
> "reduce()" member will return the original number. Given any collection of
> indices, a Domain's "is_valid()" member will return true if this
collection
> could have been obtained from a call to "expand(i)" using "i" in range
> [0,num_elements()-1], and false otherwise.
>
> For box_domain, "reduce()" reproduce the same algorithm Boost.MultiArray
> uses to locate an element in memory (using shape() and index_bases()),
> "expand()" performs the opposite task as you earlier proposed and
> "is_valid()" simply compare each element of the provided collection of
> indices to index_bases and index_bases + shape.

This is similar to the encode and decode functions in apl.  The
difference is with the direction (increasing or decreasing
w.r.t. value) of the strides and the fact that strides are already
calculated in index_list; whereas in the encode/decode functions, they
have to be calculated from scratch (from the shape vector).

Here's the equivalent to reduce and expand using functions from the
std library and variadic templates:

-{{---cut here---

enum dirs //directions(used to flag order of processing an ordered
sequence).
{ dir_fwd //forward direction.
, dir_rev //reverse direction.
};

...
      int const
    my_dir
    /**@brief
     *  Either 0 or 1.
     *    If 1, then my_strides[i]/my_strides[i+1] > 0
     *    Otherwise,  my_strides[i+1]/my_strides[i] > 0
     */
    ;
        typedef
      unsigned
    index_t
    ;
        typedef
      std::vector<index_t>
    strides_t
    ;
      strides_t
    my_strides
    /**@brief
     *  The strides for each axis in the array.
     */
    ;
...
      template
      < typename InpIter
      , typename OutIter
      >
      OutIter
    init_strides
      ( InpIter sizes_beg
      , InpIter sizes_end
      , OutIter strides
      )
      /**@brief
       *  Helper function for int_strides( dirs, Sizes...)
       */
    {
        *strides=1;
        return
        std::partial_sum
        ( sizes_beg
        , sizes_end
        , ++strides
        , std::multiplies<typename InpIter::value_type>()
        );
    }
      template
      < typename... Sizes
      >
      index_t
    init_strides
      ( dirs a_dir
      , Sizes... a_size
      )
      /**@brief
       *  Calculates strides of the array with shape, a_size...
       *  If(a_dir==dir_fwd) then strides increase with index;
       *  otherwise, they decrease with index.
       *  Returns the size of the whole array, IOW, the
       *  produce of a_size...
       */
    {
        strides_t const sizes({a_size...});
        index_t result;
        if(a_dir==dir_fwd)
        {
              auto
            it_v=init_strides
              ( sizes.begin()
              , sizes.end()
              , my_strides.begin()
              );
            result=*(--it_v);
        }
        else
        {
              auto
            it_v=init_strides
              ( sizes.rbegin()
              , sizes.rend()
              , my_strides.rbegin()
              );
            result=*(--it_v);
        }
        std::cout<<"result="<<result<<"\n";
        return result;
    }

      index_t
    rank()const
    {
        return my_strides.size()-1;
    }
...
      template
      < typename... Index
      >
      index_t
    offset_at_indices
      ( Index... a_index
      )const
      /**@brief
       *  The offset of element in my_data
       *  corresponding to indices, a_index...
       */
    {
        strides_t const indices({a_index...});
        index_t const offset
          = std::inner_product
            ( indices.begin()
            , indices.end()
            , my_strides.begin()+my_dir
            , index_t(0)
            );
        return offset;
    }

      template
      < typename InpIter
      , typename OutIter
      >
      void
    indices_at_offset
      ( index_t a_offset
      , InpIter strides_beg
      , InpIter strides_end
      , OutIter indices
      )const
      /**@brief
       *  Helper function for unary indices_at_offset(see below).
       */
    {
        std::transform
          ( strides_beg
          , strides_end
          , indices
          , [&a_offset](index_t stride)
            {
                index_t index=a_offset/stride;
                a_offset-=index*stride;
                return index;
            }
          );
    }

      std::vector<index_t>
    indices_at_offset
      ( index_t offset
      )const
      /**@brief
       *  Inverse of offset_at_indices
       */
    {
        std::vector<index_t> indices(rank());
        if(my_dir==dir_fwd)
        {
            indices_at_offset
            ( offset
            , my_strides.rbegin()+1
            , my_strides.rend()
            , indices.rbegin()
            );
        }
        else
        {
            indices_at_offset
            ( offset
            , my_strides.begin()+1
            , my_strides.end()
            , indices.begin()
            );
        }
        return indices;
    }

-}}---cut here---

In the above, indices_at_offset corresonds to index_list expand,
and offset_at_indices corresponds to index_list reduce.

>
>> Aside from runtime dimension
>> selection, what does 'box' add beyond the current extent_gen-based
>> capabilities.  (Compile-time) dimension-independent code can be
>> written by templating on the MultiArray NumDims parameter
>> (
>
http://agentzlerich.blogspot.com/2010/01/providing-fill-and-foreach-algorithms.html
>> is an example).
>
> Excluding runtime dimensionality, the sole advantages I see for data
> organized in a "box" way are:
> - it facilitate coding concerning the current position using "expand()",
> - it lighten the treatment of "boundary conditions" using "is_valid()" (By
> boundary conditions, I refer to the special treatment that is sometime
> required, in the form of "ifs", around the edges of the MultiArray.) and
> - for the cases where, although it is ok to store the data using a "box"
> organization, the same data has to be accessed in a different way that is
> not trivial to implement (The access could then be made using a different
> domain than a box_domain.).
>

I didn't notice any operations for doing transforms.  Methods for
doing transforms, as well as all the stride calculations as well as
many other operations are described in the book:

  http://web.engr.oregonstate.edu/~budd/Books/aplc/

You might find it useful.

AFAICT, one advantage that an explict dimension, which multi_array
has, over this dynamic dimension is that it allows compile-time
diagnosis of invalid number of indices.  For example, if you decided
to supply an operator[](unsigned Index), then the type of the result
would have 1 less dimension.  For example, with a multi_array<T,N>,
the result would be Multi_array<T,N-1>, if N>1, or T& if N==1.  How
would you handle this with index_list?

[snip]

-regards,
Larry