public inbox for gsl-discuss@sourceware.org
 help / color / mirror / Atom feed
* Re: containers tentative design summary
@ 2009-11-03 19:44 Gerard Jungman
  2009-11-09 20:41 ` Brian Gough
  0 siblings, 1 reply; 23+ messages in thread
From: Gerard Jungman @ 2009-11-03 19:44 UTC (permalink / raw)
  To: gsl-discuss


On Fri, 2009-10-30 at 16:03 +0000, Brian Gough wrote:
> 
> I don't think this interface is safe to use in a normal manner it it
> requires an explicit cast to pass a non-const object to a function
> accepting the const version.

You really do want to use a strongly-typed language.
But C is not it.

> Any error in the argument of an explicit
> cast is undetectable.

If the container stuff is designed properly,
then it hardly matters, because it just works
anyway, for any container type for which the
operation is meaningful.

If your world were limited to the vector and
const_vector types, as in the example, it wouldn't
matter what you did, the result would obviously be correct.

If, on the other hand, you try to jam your externally-defined
square peg into the GSL round hole, it could be bad. So?
How is this surprising?

If you think about the full container design, you will see that
it is possible to design a rational and useful system, where the
types fully inter-operate. The basis would be a generic multi-array
implementation. The vector example given was purely to address the
const-ness problem, and I think it does that.

It's not enough to scare people with notions of the wrong-cast
bogeyman. You should explain how these things can happen. When
you face the fear directly, you may realize it is not such a
big deal after all.


> We should not require the use of explicit casts
> for common operations, it would endorse a bad practice.

Actually, it endorses a "C" practice. It's "bad"-ness is
a theological issue. See the attached paper by Siff et al.,
[sorry, I tried to attach the paper, but the list won't
 accept a message with an attached pdf; the paper is
 "Coping with Type Casts in C", 1999 Bell-Labs report]
where this construct is labeled "physical sub-typing".
They give data on actual usage in the world, see Table 1.

The whole business is far from perfect. But, again, if you
want to do things in C, you need tricks, because the
language has no support for expressing relationships
between types.


In my opinion, the idea of a container library in C is
somewhat brain-damaged anyway. But here we are...


There are other ways to design containers. Some ways
involve radical change to the GSL container methodology.
For example, the const-ness of the data pointer could be
explicitly divorced from the layout meta-data which constitutes
the rest of the container. After all, the headaches with const-ness
come from the simple fact that the data pointer is embedded in the
struct and it conflates the semantics in un-desirable ways.

Example separation:

  typedef struct {
    size_t size;
    int    stride;
  } gsl_vector_layout;

  gsl_fnc_on_stuff(double * d, const gsl_vector_layout * l);
  gsl_fnc_on_const_stuff(const double * d, const gsl_vector_layout * l);


There are many ways to skin the cat. But a change like this
would effect all gsl interfaces involving container types.
It also has other noticeable imperfections.

Choose the form of the destroyer...

--
G. Jungman

^ permalink raw reply	[flat|nested] 23+ messages in thread
* Re: containers tentative design summary
@ 2009-11-20  9:37 Justin Lenzo
  0 siblings, 0 replies; 23+ messages in thread
From: Justin Lenzo @ 2009-11-20  9:37 UTC (permalink / raw)
  To: gsl-discuss

> At Thu, 29 Oct 2009 17:40:30 -0400,
> James Bergstra wrote:
> > If the purpose is to protect the user from accidentally messing around
> > with data then, as Gerard suggests, maybe we shouldn't bother.  This
> > is not a battle that we can win in C.  Good naming conventions for
> > functions, which indicate the arguments that will be modified, is the
> > most that a C library is expected to provide.
>
> The purpose is to make programs safer, rather than provide any hints
> for optimisation.  The current system is type-safe and gives a
> "discarding const" compiler warning if people try to pass const
> objects to functions as non-const arguments - these seem like useful
> features.  It's not clear to me what the actual benefit would be if we
> only had non-const vectors and matrices.
>

What if you went with a structure like the following:

typedef struct {
  const double *data;
  double *wdata;
  ...
} gsl_vector;

When a gsl_vector is allocated, the array is attached to 'wdata' and
aliased to 'data'.  When a read-only vector view is made, it returns a
new instance of the gsl_vector datatype where the 'data' and other
members are copied over, but sets 'wdata' to NULL.  Basically, you
would be allowing the same datatype to cover allocated vectors,
read-write views, and read-only views.  It seems to me this achieves
the goal of protecting underlying read-only data, as long as your
willing to check for NULL pointers in the routines that write to a
vector.  I don't know anything about the security implications, so
forgive me if this is a dangerous approach.


^ permalink raw reply	[flat|nested] 23+ messages in thread
* containers tentative design summary
@ 2009-10-05 10:12 Gerard Jungman
  2009-10-05 14:50 ` James Bergstra
  0 siblings, 1 reply; 23+ messages in thread
From: Gerard Jungman @ 2009-10-05 10:12 UTC (permalink / raw)
  To: gsl-discuss

[-- Attachment #1: Type: text/plain, Size: 184 bytes --]

A summary of the tentative design for containers. Code will follow,
as soon as a few things are sorted out. (See previous post on
questions about block/vector/matrix).

--
G. Jungman


[-- Attachment #2: gsl-containers-design-summary.txt --]
[-- Type: text/plain, Size: 11270 bytes --]

** Multi-array.

   A multi-array is a data structure managing a reference to
   a data segment (think double * as an origin) and enabling
   indexed access to the data. Indexing is controlled by
   an indexing map, which translates indices into an
   offset in the underlying data segment. The number
   of indices (dimension of the domain of the indexing map),
   is called the rank of the indexing map.

   The data segment may be externally controlled, or internally
   controlled by the multi-array. In essence, "control" means
   control of the allocation/deallocation process for the segment.

   No notion of view is needed. When the data is internally
   controlled, we can think of the multi-array as "original".
   When the data is externally controlled, we can think of
   the multi-array as a view. But this is not a precise way
   of thinking. More precisely, the life-cycle of the data
   is simply independent of the indexing semantics, and the
   need for a separate notion of view evaporates.

   This lack of a separate view type is not shared by some other
   designs. For example, boost introduces a distinct reference
   type boost::multi_array_ref, in addition to boost:multi_array.
   The reason for this choice is not clear. It may be for absolutely
   optimal efficiency. Certain choices are easier to make if one
   separates these types; however, my current feeling is that
   the extra type is not justified.

   If we decide to introduce separate view types later, we must
   insure that it is done in a way which does not force other
   interfaces in the library to grow combinatorially. And
   without doubt, the extra types must have the same generic
   interfaces as the non-view types; this is one area where
   GSL fails miserably.



** Indexing map.

   An indexing map translates an index vector to an offset. It is
   specified by a shape. At a minimum, shapes include the dimension
   of each index (number of indices addressed). When a shape allows
   an index in place i to have only one value, that index place
   is called "invariate". If it can have no values (null set), it is
   called "null". Otherwise (the normal case) it is called "variate".

   The number of indices (dimension of the domain of the indexing map),
   is called the "static rank" or simply "rank" of the map. The
   number of variate indices in a given map is called the
   "dynamic rank" or "effective rank" of the map.

   Examples using C-array notation:

     double x[3][4][5];  // static rank = 3, dynamic rank = 3

     double x[3][1][5];  // static rank = 3, dynamic rank = 2
                         // [3] and [5] are variate, [1] is invariate



** Structure of indexing maps.

   We choose the most general affine form for indexing maps.
   Therefore, the offset is calculated as
   
     offset = sum(n, index[n] * stride[n])

   The data segment is stored as a pointer to the origin of
   the data. Therefore, the most general map calculation
   takes the affine form

     location = origin + sum(n, index[n] * stride[n]).



** Slicing as composition of indexing maps.

   Let i[N] be an element of the index set for a rank N indexing map.
   Let j[M] be an element of the index set for a rank M indexing map.

   Let delta[N][M] be an integer-valued matrix and b[N] an
   integer-valued vector. Then the following formula which
   maps M-dimensional indices to N-dimensional indices
   defines a rank M index map from any rank-N index map.

     i[n] = b[n] + sum(m, delta[n][m] * j[m])

   Notice the functoriality. "Old" indices are determined in
   terms of "new" indices. In other language, the new indices
   are pulled-back to the old context. This is the correct
   functoriality to allow arbitrary chained composition of
   index maps, and therefore slices of slices, to any order.

   It should be clear that the index map generated this way
   is always affine, and all affine indexing maps can be
   generated in this way.

   Therefore, a general slice of an affine indexing map
   is specified by a vector b[] and matrix delta[][], of
   the indicated dimensions.

   The reader can convince himself that any form of strided
   access is supported in this way. For example, it is equally
   easy to slice out the super-diagonal of a matrix as it is to
   slice out a row or column, or for that matter to slice out
   any regularly strided sub-matrix, etc.

   The underlying direct data access, by computation of an
   offset, is itself an affine indexing map, specialized to
   the case where the "old" index is one-dimensional. So
   the design space of affine indexing maps is closed.



** Multi-array implementation requirements.

 * Indexing calculations must be as fast as possible. The indexing
   arithmetic should be inlined when possible. The shape data
   structure should enable the fastest possible computation,
   consistent with general strided access.

   This suggests that different rank objects be of different static
   types, since this is the most straight-forward way to insure that
   the correct indexing calculation is chosen at compile-time. This
   design corresponds to the usual designs (such as in GSL), where
   vectors and matrices (for example) are different static types.
   This choice also corresponds to the design of boost::multi_array,
   where the rank is a template parameter and therefore a
   compile-time constant which generates a separate type.

 * The shape of a multi-array should be dynamic. This corresponds
   to the usual design, such as in GSL, where the dimensions of 
   a matrix (for example) are run-time constants, not compile-time
   constants.

 * It should be possible to construct multi-arrays without
   heap allocations. For example:

     double data[64];
     gsl_marray_2 ma_8x8 = gsl_marray_2_build(data, 8, 8);

   Placing such objects on the stack has good implications
   for efficiency and is very natural from the client's point
   of view. For example, in C++, one simply writes

    std::vector<double> v;

   which places the object on the stack as an automatic variable.



** Multi-array implementation choices. These are open for debate,
   but these are the choices I have made. No loss of functional
   generality is implied by these choices, only some loss of
   syntactic generality.

 * Multi-array indexing conforms precisely to C-style array
   indexing. Shapes are completely specified by dimensions and
   strides. This means that indices always "start" at 0, and
   end at dimension-1 (inclusive).

 * The notion of storage order does not enter. The "fast" index
   is always to the right, conforming to C array semantics. See
   the discussion of matrices below.



** Vector.

   A vector manages a data segment and supports indexed access
   to the data through a rank one indexing map. However,
   vectors may be required to satisfy extra constraints, not
   compatible with rank one multi-arrays. Also, they may provide
   other capability, such as indexing which "starts" at some
   non-zero value. Personally, I have no desire for this extra
   sort of functionality. It would only add overhead to an
   otherwise clean design. My first inlination is to
   implement vectors as follows:

       typedef  gsl_marray_1  gsl_vector; // DONE

   However, the reality may not be quite as simple as this in
   the end. However it is done, I think it is a good idea
   to preserve the name "vector" in the design, although
   it will probably imply some significant boiler-plate
   in function declarations, since C does not allow us
   to simply inherit the member function namespace
   of another class. SIGH.



** Matrix.

   Matrices are clearly different from rank two multi-arrays.
   There are necessary extra capabilities, and there are also
   constraints, in order to make them conform to the notions
   used in linear algebra packages.

   The extra capability corresponds to the identification of "row"
   and "column" indices. Multi-array semantics are governed by data
   locality, so the notion of "fast" index is immutable. But a matrix
   may either identify its rows with fast index traversal (C order)
   or its columns with fast index traversal (Fortran order). Matrices
   will need to remember which and to adjust their indexing computations
   accordingly.

   The extra constraint on matrices, in comparison to multi-arrays,
   is the need for the fast index to traverse the data with stride 1.
   For example, BLAS interfaces allow for vectors of arbitrary
   stride, but do not allow for a matrix "fast index" stride.

   Example:
     Let A be a 3x3 matrix. Let B1 be a 2x2 matrix defined as a
     view (sub-matrix) of A defined by the x's in the following
     diagram:
     
       x x -
       - - -
       x x -

     Then B1 is a matrix, in the BLAS sense. It can be passed
     to BLAS functions without further reshaping. To be exact,
     it would be passed with origin equal to the given origin
     and with tda = 6 and dimensions (2,2).

     Let B2 be the 2x2 sub-matrix view defined by the following
     diagram.
     
      x - x
      - - -
      x - x

     Then B2 is _not_ a BLAS compatible matrix, since there is no
     parameter in the BLAS interfaces which allows for the non-unit
     stride in the fast index. Some reshaping, involving data
     motion/copy, would be required to pass this view to BLAS
     functions.

   Therefore, I do not think that generic rank two multi-arrays can
   be identified with matrices. Rather, a lightweight conversion
   from rank-two multi-arrays to matrices must be provided, and this
   conversion must fail at run-time if the multi-array fast index
   stride is not unity.



** Tensor

  No explicit notion of tensor currently exists, beyond that
  expressed in the tensor/ extension package for GSL, by
  Jordi Burguet-Castell.

  The Burguet-Castell package does not meet the requirements
  layed out above, mainly because it looks and works too much
  like the GSL types, having been modelled on them for obvious
  reasons. In particular, it allows the rank to be determined
  at run-time. This may or may not be desirable, depending on
  how we want tensors to interact withy multi-arrays.

  However, we draw some motivation from the Burguet-Castell package.
  A tensor is defined to be a multi-array where the shape is constrained
  to have all dimensions equal.

  Like in the case of matrices, we encounter some small technical
  problems in the identification with multi-arrays. Since the
  shape of a multi-array is dynamic, the compatibility of
  a given multi-array with the tensor concept is question
  that must be deferred to run-time.

  The question of the role of multi-arrays in the implementation
  of tensors remains open, as it does for matrices. But we do
  require a certain level of interoperability, perhaps expressed
  as light-weight conversion functions which fail to convert a
  multi-array to a tensor if the shape is not consistent. This
  seems like the most likely solution. It also allows the tensor
  implementation to be independent of the high-level multi-array
  interfaces, which strikes me as the right way to do it, from
  both a general design point of view and an efficiency point of
  view.


^ permalink raw reply	[flat|nested] 23+ messages in thread
end of thread, other threads:[~2009-11-20  9:37 UTC | newest]

Thread overview: 23+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2009-11-03 19:44 containers tentative design summary Gerard Jungman
2009-11-09 20:41 ` Brian Gough
2009-11-09 23:06   ` Gerard Jungman
2009-11-14 15:25     ` Brian Gough
2009-11-15  9:13       ` Tuomo Keskitalo
2009-11-15 16:44         ` Jonathan Underwood
2009-11-15 18:41           ` Robert G. Brown
2009-11-16 11:56         ` Brian Gough
  -- strict thread matches above, loose matches on Subject: below --
2009-11-20  9:37 Justin Lenzo
2009-10-05 10:12 Gerard Jungman
2009-10-05 14:50 ` James Bergstra
2009-10-05 23:00   ` Gerard Jungman
2009-10-05 23:45     ` James Bergstra
2009-10-06 19:59       ` Gerard Jungman
     [not found]     ` <645d17210910060537s762d6323pfd2bec8590ad28e9@mail.gmail.com>
2009-10-06 20:02       ` Gerard Jungman
2009-10-23 21:28     ` Brian Gough
2009-10-27 23:06       ` Gerard Jungman
     [not found]         ` <7f1eaee30910271628h70785125m68e47c7a7b5c25b7@mail.gmail.com>
2009-10-27 23:49           ` Gerard Jungman
2009-10-29 18:06         ` Brian Gough
2009-10-29 20:41           ` Gerard Jungman
2009-10-29 21:40             ` James Bergstra
2009-10-30 16:54               ` Brian Gough
2009-10-30 16:54             ` Brian Gough

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).