Le mercredi 21 mai 2014 14:55:34 UTC+2, Tim Holy a écrit :
>
> I agree that vectorizing is very nice in some circumstances, and it will
> be
> great if Julia can get to the point of not introducing so many
> temporaries.
> (Like Tobi, I am not convinced this is an easy problem to solve.)
>
> But there are plenty of cases where loops are clearer even when
> vectorization
> is both possible and efficient. For example, surely you can't tell me that
>
> L = similar(A)
> m,n = size(A)
> L[2:m-1,2:n-1] = A[3:m,2:n-1] + A[1:m-2,2:n-1] + A[2:m-1,3:n] +
> A[2:m-1,1:n-2] - 4A[2:m-1, 2:n-1]
>
> is easier to read and understand than
>
> L = similar(A)
> for j=2:size(A,2)-1, i=2:size(A,1)-1
> L[i,j] = A[i+1,j] + A[i-1,j] + A[i,j+1] + A[i,j-1] - 4A[i,j]
> end
>
> I completely agree. When I was talking about vectorized code, I meant what
Andreas showed.
> Not only does the loops-version have fewer keystrokes, it's a heck of a
> lot
> clearer that what you're doing is computing the 2d discrete laplacian.
>
> --Tim
>
> On Wednesday, May 21, 2014 05:09:57 AM Andreas Lobinger wrote:
> > Hello colleague,
> >
> > On Wednesday, May 21, 2014 7:33:01 AM UTC+2, Andreas Noack Jensen wrote:
> > > Please consider b += Diagonal(c) instead of diagm. Diagonal(c) only
> stores
> > > the diagonal elements but works like diagm(c) for matrix arithmetic
> > >
> > > Vectorized code is always easier to understand
> > >
> > >
> > > That is a strong statement. I have vectorised MATLAB code with repmat
> and
> > > meshgrids that is completely unreadable, but would be fairly easy to
> > > follow
> > > if written as loops. I really enjoy that I can just write a loop in
> Julia
> > > without slow execution.
>
Me too :) I agree about my statement: I completely overlooked complicated
slicing and stuff. But I also said that a part of what made explicit loops
less easier to understand are explicit indices. Same thing goes for
"vectorized" code.
My whole point is to say that, in general, one should divide the problem
into understandable pieces to write clear code. Vectorized expression are a
way to get that. For instance,
A -= Diagonal(B)
is very easy to understand, if not obvious. You almost wrote in plain
english there, and that's a strength: no need to *actually* understand the
code to grasp what it does.
for i=1:length(B)
A[i, i] -= B[i]
end
is not that hard to understand because it's really short, but there, the
only way you can get the intent of the programmer is by browsing the whole
code and understanding the whole thing.
To get back to Tim example.
On could think of creating a DiscreteLaplacian2D function taking a matrix A
as an input. This is actually what I would do to keep the complexity of any
given piece of code very low. Let's solve the heat equation iteratively,
for that we need to calculate:
DiscreteDerivative(A, 1) - alpha.*DiscreteLaplacian2D(A, [2,3])
or even better, knowing the shape of A:
TimeDerivative(A) - alpha.*SpaceLaplacian2D(A)
are very easy to grasp and understand. I've not directly used loops : this
is plain vectorized code. It's obvious and easy to maintain.
I'm sure you will agree. But what's the problem then? This is not *really*what
you meant about vectorized code.
Well, let's handle the same problem with the current trend : this is not as
fast as Julia could, this might be fast enough but definitely not the
fastest. It's also using more RAM that it should, so why not *rewriting it
into loops*? Yes, of course, why not! The only problem is that the only way
to achieve something usefull would be to solve the *whole* equation in one
fat bunch of loops. Which means:
- Duplicated code all over the place. You'll have to write/(copy/paste)
your devectorized laplacian for *each* single equation that needs it.
- Unreadable code.
- Longer code.
- More places to make mistakes.
- More places to search for to correct copy/paste mistakes.
- Bigger mental map for you project.
- Fewer external contributors.
- more "end"s ;)
>
> > If you write code in Matlab (or similar) that's using repmat (nowadays
> > bsxfun) and meshgrids (and reshape) you are operating on wrong side of
> > vectorising your code. Most likely (and there might be good reasons for
> > that) your data is not in the 'right' shape, so you reorder more than
> you
> > calculate. For problems like this explicit loops are often the better
> way
> > (so straightforward in julia) OR you take some thinking, why the data
> needs
> > to be reordered for certain operations.
> >
> > The right side of vectorized code is the SIMD thinking which does
> > operations (some times linalg) on Matrix/Vector forms, and that in the
> best
> > sense leads to one-liners like
> >
> > M += ((A*B) + (d*C) + (E .* M0))
> >
> > which are easy to read/write/understand.
> >
> > Expressions like this seems to create in julia an incredible ammount of
> > intermediates and temporary allocations, while actually the output array
> > (M) is already allocated and the + and * operations just needed to be
> > executed with the 'right' indexing.
> >
> > Some people tend to solve this by writing explicit loops and create half
> a
> > page of code for the above one-liner.
>
