Re: Octave 3.0

[Etienne Grossmann]

  Hello,

  this is just to add my 2c to the multi-dimensional array
discussion. In fact, these 2c come from the manpage of pdl (Perl Data
Language - pdl.sourceforge.org), which support arbitray number of
dimensions, slices that are not copied, dimension reshuffling etc.

======================================================================
NAME
       PDL::Slices -- Stupid index tricks

SYNOPSIS
         use PDL;
         $a = ones(3,3);
         $b = $a->slice('-1:0,(1)');
         $c = $a->dummy(2);


DESCRIPTION
       This package provides many of the powerful PerlDL core index manipula­
       tion routines. These routines are usually two-way so you can get a
       unit matrix by

        $a = zeroes(1000,1000);
        $a->diagonal(0,1) ++;

       which is usually fairly efficient. See PDL::Indexing and PDL::Tips for
       more examples.

       These functions are usually two-way:

        $b = $a->slice("1:3");
        $b += 5;               # $a is changed!

       If you want to force a copy and no "flow" backwards, you need

        $b = $a->slice("1:3")->copy;
        $b += 5;               # $a is not changed.
======================================================================

  Cheers,

  Etienne

On Mon, Mar 03, 2003 at 09:49:32AM -0500, Lippert, Ross A. wrote:
# (this post is long than I thought it would be.
# contents:
#   re: zero-copy reshape (and its myriad uses)
#   re: sparse matrices
#   re: matlab compatability
#   re: MPI and PVM
# )
# 
# re: zero-copy reshape
# >reshape as currently implemented requires a matrix copy,
# >but that would be easy to fix.
# I agree that a copy is inconvenient.  Usually it is not an issue
# since the other ops are more expensive than the copy.
# 
# But, one thing that I think could be really cool on this topic
# is to abstract the block of doubles allocated to a matrix/vector
# from the way the block is indexed, supporting multiple
# views for the same data, e.g.:
#  B = A'
# or
#  B = reshape(A,prod(size(A)),1)
# would allocate no new memory for B, just create a new 'view'.  A new
# block only gets allocated when one writes to an entry of B.  This
# kind of view business should be compatable with the basic BLAS ops
# since they allow transpose bits and strides to be set appropriately.
# 
# I noticed that GSL supports this trick, to the extent that
#  v = A(1,:)
# is just another view.  I'm not sure if they handle
#  B = A([1 3 4],:)
# as a view, and I don't see how one could do a gemm on such a B
# without forcing a tmp copy.
# 
# BTW I am not an advocate of GSL-ing octave.  I've looked at some
# of the GSL sources, and I'm not as impressed with the implementation
# as I am with the interface.
# 
# >There is something to be said for convenient data
# >structures, especially if it is convenient to operate
# >with slices.  The results are going to be more readable
# IMHO I'd change 'especially' to 'only'.  Sometimes the operations
# we think we want aren't the right ones, and when we figure out what
# ops we want, we figure out better representations for our data.
# 
# >than tricky indexing and reshaping operations.  It's
# >not clear to me that Matlab does the right thing --- not
# >while squeeze/reshape are required to transform a
# >slice using a matrix operation.
# Definitely.  I think it is very hard to work out what one
# really should have basic ops on multidimensional data.  One
# ugly, but sensible, thing one should have is something like
#  C = CONTRACT(A,B,[2 3],[4 1])
# where C(a,b,c) = SUM_{xy} A(a,x,y)*B(y,b,c,x)
# but that's kind of nasty.  On the other hand one also can do
#  C = reshape( reshape(A,na,nx*ny)*reshape( reshape(B,ny*nb*nc,nx)' 
,nx*ny,nb*nc) ,na,nb,nc)
# which is also kind of nasty, but at least it introduces no new ops
# to the language.
# 
# And yes, you have to have the cock-eyed view of 
# matrices like I have , so you are right that this is for the
# sophisticated user and not for everyone.  Maybe being forced
# to do things like this makes one sophisticated.  Maybe I'd just
# too proud of myself for tricks like this.
# 
# If zero-copy reshape/transpose were implemented, and some support
# for inline functions were around, then perhaps one could just
# write .m scripts to support N-D arrays and sparse matrices.
# 
# re: sparsity
# For sparse matrices, I agree that idxop is definitely for sophisticates
# and wrapping it up might require making some sort of special 'struct'
# with some kind of overloading.  I'm kind of partial to some special
# OO struct type which has a special field which says 'I have functions
# in me', and then when someone does A*B and the interpreter sees that one
# of them is a struct, and instead of the usual 'I don't know what to do with
# structs' error message, it can check for the OO flag and look in the struct
# fields for 'times' to pull out the appropriate functions.  This is a
# nice way to take advantage of the fact that functions are first-class
# (or nearly first class) objects in octave, and I think it is better than
# MATLABs "greedy overloading by directory name" crap.
# 
# function s = sparse(A)
#   [i,j,a] = find(A);
#   s.OO = true;
#   s.i = i;  s.j = j; s.a = a;
#   s.times = sparse_times;
#   s.plus  = sparse_plus;
#   ...
# endfunction
# 
# BTW just one last thing on idxop.  The thing that really necessitates it
# is the fact that += doesn't do 'the right thing' (which may not be right
# to some people) when the LHS is multiply indexed:
#  v = zeros(2,1);
#  v([1 2 2]) += 5
# currently sets v = [5 5] instead of [5 10].  A sparse multiply could then
# be
#  [i,j,a] = find(A);
#  y = zeros(size(A,1),1)
#  y(i) += a.*x(j)
# I had this argument before with some people on this list (and I'm content tha
# my position may not be the right one), and I am aware that it would require
# some serious fiddling with octave sources to do it my way.  Then there is
# the issue of not having a max= or a min= operator for the other idxops,
# though one might co-opt &= and |= for these services.
# 
# re: MATLAB compatability
# The bit about rewriting code on the web is a really good one, for
# which I don't have a good answer.  MATLAB has allowed people to
# write some obnoxious things, and for octave to support that code
# it needs to support some of these obnoxious things.  Perhaps, it
# can be done with little change to octave other than some mods to
# the basic operations then wrapped up in .m's.  That would be nice.
# 
# My point is that octave presents an opportunity not just to be
# compatable with MATLAB, but to go that extra mile and do things
# smarter than MATLAB.  (BTW another personal example of this is that
# octave includes the LAPACK sylvester equation solver, which, to my
# shock, MATLAB doesn't, and I was solving sylvester equations a lot
# at one point and probably will need to again in the future).  I'd
# like to see new features put in because it is a good idea, whether
# from MATLAB or elsewhere, with a MATLAB-like interface when
# applicable.
# 
[snip]


-- 
Etienne Grossmann ------ http://www.isr.ist.utl.pt/~etienne