On Thu, Jun 25, 2009 at 10:55 AM, John W. Eaton
<address@hidden> wrote:
On 25-Jun-2009, Judd Storrs wrote:
| What about a new type instead? In this form it seems like a small project I
| can work on by following diag as an example (I have limited time for this,
| but Robert Short's example is inspiring). I'm having trouble thinking of
| *the* awesome name for this and the name isn't really important for now, but
| what about something like adapt(), expand(), grow(), match(), merge() or
| join()? The synyax would look like this:
|
| X / adapt(sum(X,1))
| adapt(sum(X,2)) \ X
|
| Would something like this be acceptable or useful to anyone except me?
Instead of a new type, what is wrong with just defining a function to
do the special operation?
Yes, you could write a bazillion specialty functions.
Maybe this is a difference in perspective based on the type of data I process. If you're thinking about things in terms of linear algebra and 2D matrix operations, diag() is the correct abstraction for scaling rows and columns. For matrix operations the advantages of diag() beyond performance are that it is mathematically expressive in the sense that a matrix _expression_ maps pretty much verbatim into octave syntax. When you're working with 2D matrices you're generally thinking about matrix multiplication and matrix inversion and diag() fits right in. In these cases you end up with an octave _expression_ that *looks* like the math.
Things are different when you approach working with 3D and 4D volumes, IMHO. Most of the stuff I do is thinking working on slices and planes. i.e. translating tensor type expressions:
X_ijk = (A_ij * B_k + C_ik) / (D_j + sum_all_j E_jk)
Of course the natural solution is to wrap this _expression_ into for loops. As you know, this looks like the math but performs terribly:
X(i,j,k) = ( A(i,j) .* B(k) + C(i,k) ) .* ( D(j) / sum(E(:,k)) )
You can eat up memory and expand all of A,B,C, and D to full size. Then
sumE = sum(E,1)
X = (A.*B + C).*(D ./ sumE(k))
But this can eat memory very fast for large arrays. You can instead set up your dimensions correctly so that you can use bsxfun:
# A is Nx by Ny by 1
# B is 1 by 1 by Nz
# C is Nx by 1 by Nz
# D is 1 by Ny by 1
# E is 1 by Ny by Nz
X = bsxfun(@rdivide,bsxfun(@times,bsxfun(@plus,bsxfun(@times,A,B),C),D),sum(E,2))
But this is beginning to be unreadable and often there are faster ways. To get something that actually performs very well in octave, you end up doing all sorts of reshape, transpose, reshape, bsxfun, permute gymnastics and the resulting code doesn't really look like the math at all. You could put any solution in a function or write an .oct, and then use
X = compute_equation_42(A,B,C,D,E)
But that doesn't really preserve the mathematical _expression_ either.
With the proposed type and the matrices dimensioned as for bxsfun the following could be expressive and efficient
X = ( adapt(A) .* adapt(B) + adapt(C) ) .* adapt(D) ./ adapt(sum(E,2))
Or (more likely) if A, B, C, and D were created as adapt type from the start then
X = (A .* B + C) .* ( D ./ sum(E,2) )
Also, I'm not a convinced as you are that modifying the element-wise builtins to have bsxfun behavior would be a problem. The only cases where the enhanced versions would succeed is when the dimensions of the arguements meshed in a bsxfun compatible way. It turns out that many octave/matlab functions already generate bsxfun compatible results (i.e. sum, mean, std, etc.).
--judd