octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

octave "Most Wanted" feature


From: robert bristow-johnson
Subject: octave "Most Wanted" feature
Date: Mon, 27 Nov 2006 09:56:06 -0500

i am resending this to the octave-maintainers as per Dr. Eaton's request.

i know it's long, i know it has been written off by Cleve Moler and the 
MathWorks. but i am hoping that Octave doesn't repeat the same mistake.

this really is important, and DOABLE in an open source resource like Octave, 
but even though i can code pretty well in C and can dabble in C++, i cannot 
crack the whole source code of the project.  perhaps someone more knowledgable 
on how the project is built can.  here is my letter:


Dear Dr. Eaton,

I know this is six years late (since my diatribe at comp.soft-sys.matlab
that shows you also in the thread), but I thought I might again dredge
this issue up with regard to Octave in lieu of MATLAB. Since I can't
edit the Octave Wiki and I do not see exactly how to get on the octave
devlopers mailing list, I thought I would send this to you regarding
what I think is an important improvement to Octave that would make it
significantly better than MATLAB for signal processing engineers.  (I'll
also confess that I am now starting to use a linux machine and Octave
has become my math tool replacing MATLAB.)

This is, if you recall, about the array base index issue that I have
brought up a few times at comp.soft-sys.matlab (with cross-posting at
comp.dsp).

Now i realize that Octave is not MATLAB, but I would like to repeat my
basic indictment against this hard-wired 1-based indexing in MATLAB here
for the hope that Octave need not continue the mistake (you might want
to print this out - and use a monospaced font):

______________

 From Getting Started with MATLAB, v. 5:

"MATLAB is a high-performance language for technical computing.  It
integrates computation, visualization, and programming in an easy-to-use
environment where problems and solutions are expressed in familiar
mathematical notation."

Or from an older MATLAB reference guide (ca. 1992):
"MATLAB integrates numerical analysis, matrix computation, signal
processing, and graphics in an easy-to-use environment where problems
and solutions are expressed just as they are written mathematically - ... "

Note the phrases in claims: "familiar mathematical notation" and "just
as they are written mathematically".  I submit that those claims are
false in a sense that is salient particularly for those of use who use
MATLAB for (digital) signal processing.

In Discrete-Time Signal Processing (otherwise known as "Oppenheim and
Schafer") p. 531, Eqs. (8.57) and (8.58), the Discrete Fourier Transform
(DFT) and inverse DFT are defined as:

                  N-1
     X[k] =       SUM{ x[n] * exp(-j*2*pi*n*k/N) }
                  n=0

                  N-1
     x[n] = (1/N)*SUM{ X[k] * exp(+j*2*pi*n*k/N) }
                  k=0

Whereas MathWorks defines these for MATLAB to be:

                   N
     X[k] =       SUM{ x[n] * exp(-j*2*pi*(n-1)*(k-1)/N) }
                  n=1

                   N
     x[n] = (1/N)*SUM{ X[k] * exp(+j*2*pi*(n-1)*(k-1)/N) }
                  k=1

Likewise on p. 21 Eq. (2.39) O&S define discrete convolution as:

              +inf
     y[n] =   SUM{ x[k] * h[n-k] }
             k=-inf

which is equivalent to:

              +inf
     y[n] =   SUM{ x[n-k] * h[k] }
             k=-inf

If the impulse response were causal and finite (i.e. a causal FIR where
N =# of taps), that equation would degenerate to:

              N-1
     y[n] =   SUM{ x[n-k] * h[k] }
              k=0

But MATLAB says it's this:

               N
     y[n] =   SUM{ x[n-1-k] * h[k] }
              k=1
or
              N-1
   y[n+1] =   SUM{ x[n-k] * h[k+1] }
              k=0

and, in fact, makes excuses for that in the Signal Processing Toolbox
manual (ca. 1993) p. 1-12:
"Note: The filter coefficients start with subscript 1 rather than 0.
This reflects MATLAB's standard indexing scheme for vectors."
and p. 2-56:
"The series is indexed from n+1 and k+1 instead of the usual n and k
because MATLAB vectors run from 1 to n instead of 0 to n-1."

The DFT and convolution summations are not the only examples (the way
IIR coefficients are ordered is another that occurs to me), but are *so*
fundamental that for a company whose product expresses mathematical
instructions "just as they are written mathematically", lot's of red
flags should have been flying up when they were essentially being forced
to change or contrive those well established expressions into their
fixed 1-based indexing scheme.  The fact they didn't do this 10 years
ago is no excuse for not fixing the problem now ASAP, because their
product has already become the defacto standard used by the signal
processing community to concisely exchange ideas and algorithms (in
addition to modeling and emulating these algorithms for themselves).

______________

And here is a comment from another DSP heavyweight, Eric Jacobsen, the
"Minister of Algorithms" at Intel:

> Legacy of a bad idea, in my opinion, is a poor excuse for continuing
> along ad nauseum with the bad idea.  At some point it is a benefit to
> correct the bad idea, and the longer one waits to do it the more
> painful it is.  Most of us have probably dealt with this in our
> engineering careers in some way or other, and it is rather irritating
> to see an otherwise premier outfit and product like Matlab and the
> Mathworks refuse to address it.

To which I responded:

> ...I might note is that I am proposing an extension THAT PRESERVES
> LEGACY.  These c.s-s.m guys can't squawk about that.  In addition,
> I don't want to change the base for *all* arrays with some global
> SINDEX.  I want each and every MATLAB array to have its own user-
> definable 'SINDEX' (I would call it something else) for each
> dimension of the array and that has default value of 1 when the
> array (or "matrix") is first created.  That will preserve backward
> compatibility and allow for the generalization we need.

______________

After finally getting a response from Cleve Moler, I posted what this
backward-compatible extension would mean for some matrix operations:

Not being a MathWorks insider (I can't imagine why not :-) I have to
guess a little at the structure of a MATLAB variable:

enum MATLAB_class {text, real, complex};
// I don't wanna cloud the issue considering other classes.

typedef struct
  {
  void* data;     // pointer to actual array data
  char* name;     // pointer to the variable's name
  enum MATLAB_class type; // class of MATLAB variable (real, complex,...)
  int num_dimensions;     // number of array dimensions >= 2
  long* size; // points to a vector with the number of rows, columns,etc.
  } MATLAB_variable;

name[32];       // MATLAB names are unique to 31 chars
size[num_dimensions];

if (type == text)
  {
  char data[size[0]*size[1]*...*size[num_dimensions-1]];
  }
  else if (type == real)
   {
   double data[size[0]*size[1]*...*size[num_dimensions-1];
   }
  else if (type == complex)
   {
   double data[2][size[0]*size[1]*...*size[num_dimensions-1]];
   }

When an element, A(n,k), of a 2 dimensional MATLAB array A is accessed,
first n and k are confirmed to be integer value (not a problem in C),
then confirmed to be positive and less than or equal to size[0] and
size[1], respectively.  It those constraints are satisfied, the value of
that element is accessed as:

data[(n-1)*size[0] + (k-1)];

For a 3 dimensional array, A(m,n,k), it would be the same but now:

data[((m-1)*size[1] + (n-1))*size[0] + (k-1)];

//// Note to John: I realize that Octave reverses the order of which
//// subscript varies fastest.  Octaves left-most subscript varies
//// fastest whereas MATLAB's (and C's) right-most varies fastest.
//// That's just a convention, which does not affect this issue.

I realize that the pointer to "data" can be judiciously offset so that
the subtraction of 1 from the MATLAB indices to create the C indices
would not be necessary.  I think any modern Fortran does this.  But that
is the case for *only* 1xN sized vectors.  When it's a 2-dim or 3-dim
array, this "1" base index MUST be subtracted from all of the subscripts
before the net index to this ultimate 1-dimensional C array is accessed.
  Also I realize that the MATLAB variable structure may have other
internal fields that are not described above and that I don't know
about, but I don't see any reason what that would affect the issues here.

What is proposed is to first add a new member to the MATLAB variable
structure called "base" which is a vector of the very same length
(num_dimensions) as the "size" vector.  The default value for all
elements of the base[] vector would be 1 with only the exceptions
outlined below. This is what makes this backwards compatible, in the
strictest sense of the term.

typedef struct
  {
  void* data;
  char* name;
  enum MATLAB_class type;
  int num_dimensions;
  long* size;
  long* base; // points to a vector with index base for each dimension
  } MATLAB_variable;

name[32];
size[num_dimensions];
base[num_dimensions];

Now immediately before each index is checked against the bounds for that
dimension ( >= 1 and <= size[dim] where 0<dim<num_dimensions), the base
for that particular dimension (base[dim]) is subtracted from the index
(instead of subtracting the constant 1) and then the bounds comparison
is made, and the element is accessed.

So now the value of an element of A(n,k) is accessed as:

data[(n-base[1])*size[0] + (k-base[0])];

For a 3 dimensional array, A(m,n,k), it would be the same but now:

data[((m-base[2])*size[1] + (n-base[1]))*size[0] + (k-base[0])];

Since the default is 1, this will have no effect, save for the teeny
amount of processing time need to look up the base index over using the
constant base index of 1, on existing MATLAB legacy code.

Okay, how someone like myself would use this to do something different
is that there would be at least two new MATLAB facilities similar to
size() and reshape() that I might call "base()" and "rebase()",
respectively.  Just like MATLAB size() function returns the contents of
the size[] vector, base() would return in MATLAB format, the contents of
the base[] vector.  And just like reshape() changes (under proper
conditions) the contents of the size[] vector, rebase() would change the
contents of the base[] vector.  Since rebase() does not exist in legacy
MATLAB code (oh, i suppose someone could have created a function named
that, but that's a naming problem that need not be considered here),
then there is no way for existing MATLAB programs to change the base
indices from their default values of 1 making this fix perfectly
backward compatible.

Now just as there are dimension compatibility rules that exist now for
MATLAB operations, there would be a few natural rules that would be
added so that "rebased" MATLAB arrays could have operations applied to
them in a sensible way.

ARRAY ADDITION and SUBTRACTION and element-by-element ARRAY
MULTIPLICATION, DIVISION, POWER, and ELEMENTARY FUNCTIONS:

Currently MATLAB insists that the number of dimensions are equal and the
size of each dimension are equal (that is the same "shape") before
adding or subtracting matrices or arrays.  The *one* exception to that
is adding a scaler to an array in which a hypothetical array of equal
size and shape with all elements equal to the scaler is added to the
array.  The resulting array has the same size and shape as the input arrays.

The proposed system would, of course, continue this constraint and add a
new constraint in that index bases for each dimension (the base[]
vector) would have to be equal for two arrays to be added.  The
resulting array would have the same shape and base[] vector as the input
arrays.

"MATRIX" MULTIPLICATION:

A = B*C;

Currently MATLAB appropriately insists that the number of columns of B
are equal to the number of rows of C (we shall call that number K).  The
resulting array has the number of rows of B and the number of columns of
C.  The value of a particular element of A would be:

           K
A(m,n) = SUM{ B(m,k) * C(k,n) }
          k=1

The proposed system would, of course, continue this constraint and add a
new constraint in that index bases must be equal for the dimension where
the lengths must be equal.  That is the number of columns of B are equal
to the number of rows of C and the base index of the columns of B are
equal to the base index of the rows of C.  The resulting array has the
number of rows of B and the number of columns of C and  the base index
of the rows of B and the base index of the columns of C.  The value of a
particular element of A would be:

         base+K
A(m,n) = SUM{ B(m,k) * C(k,n) }
         k=base

where base = base[0] for the B array and base[1] for the C array which
must be the same number.

Both of these definitions are degenerations of the more general case where:

          +inf
A(m,n) = SUM{ B(m,k) * C(k,n) }
         k=-inf

where here you consider B and C to be zero-extended to infinity in all
four directions (up, down, left, and right).  It's just that the zero
element pairs do not have to be multiplied and summed.

Probably matrix powers and exponentials (on square matrices) can be
defined to be consistent with this extension of the matrix multiply, but
I can deal with it at the moment.

CONCATINATION:

This would also be a simple and straight-forward extension of how MATLAB
presently concatinates arrays.  When we say:

A = [B C];

The number of rows of B and C must be equal, but the number of columns
of B and C can be anything.  The first columns of A are identical with
the columns of B and then also must the indices of those columns.  And
independent of what the column indices of C are, they just pick up where
the column index of B left off.  This rule extension defaults to what
MATLAB presently does if B and C are both 1-based arrays.  A similar
rule extension can be made for A = [B ; C];  In all cases the upper left
corner of A is identical to the upper left corner of B, both in value
but also in subscripts (so A(1,1) becomes B(1,1) just like it does now
in MATLAB).

MATRIX DIVISION ('/' and '\'):
I have to think about that a little, but I'm pretty sure a backward
compatible extension to the operations can be figgered out.  If not, it
would be an illegal operation unless the bases were 1.

FUNCTIONS THAT RETURN INDICES (min(), max(), find(), sort(), ind2sub(),
and any others that I don't know about):
It must be internally consistent (and certainly can be made to be).  The
indices returned would be exactly like the 1-based indices returned
presently in MATLAB except that the base index for the corresponding
dimension (that defaults to one) would be added to each index.  That is,
just like now in MATLAB, if:

[max_value, max_index] = max(A);

This must mean that A(max_index) is equal to max_value.

I think that this is easy enough to define.  The only hard part is to
identify all MATLAB functions that search through an array and modify
them to start and end at indices that might be different than 1 and
size[dim] as are the search bounds today.  It would instead search from
base[dim] to size[dim]+base[dim]-1 which would default to the current
operation if the base index equals one.

FOR ALL OTHER MATLAB OPERATIONS, until a reasonable extended definition
for non 1-based arrays is thunk up, MATLAB could either bomb out with an
illegal operation error if the base is not 1 or could, perhaps, ignore
the base. Either way it's still backwards compatible.

______________


Finally, Dr. Eaton, what I am asking for is your consideration of this
generalization of indexing in Octave that would make it *better* than
MATLAB.  Eventually, I predict that DSP engineers would really catch on
to this generalization (it is *not* just generalizing to 0-based arrays,
but arrays that can have *any* base index, including negative indices)
and in a decade might see such widespread use that people will wonder
"why wasn't this done earlier?".  As Eric Jacobsen said:

     "Legacy of a bad idea, in my opinion, is a poor excuse for
      continuing along ad nauseum with the bad idea.  At some point
      it is a benefit to correct the bad idea, and the longer one
      waits to do it the more painful it is."

This rigid 1-based indexing is NOT a good idea.  We need the ability to
define arrays with whatever real integer base we want for any dimension
in the array.  Please consider supporting this backward-compatible
enhancement to Octave.

Being an electrical engineer (not a comp-sci graduate), I found the
total source to Octave to be so formidable, I would not know where to
even look, but if someone more familiar with the project could spend a
little time showing me where, in the source, some of this essential
index manipulation is done and how to make the whole project, I have GCC
and would be willing to put in some coding effort.  But I am not capable
of just downloading the source (I *did* do that) and cutting right into
the patient.  Can this be made into a topic or issue in the Octave wiki
or discussion groups (wherever they are)?

Thank you for reading this, I know it is long (I hope you printed it
out).  I can also be contacted at address@hidden .  Feel free
to ask any questions.

Best regards,

r b-j

--

Robert Bristow-Johnson         address@hidden
Kurzweil Music Systems     781/890-2929 x244
Young Chang Research & Development Institute
1432 Main St.    Waltham    MA    02541-1623


--

r b-j                  address@hidden

"Imagination is more important than knowledge."




reply via email to

[Prev in Thread] Current Thread [Next in Thread]