Stefan, sum(foo) did not work for me. It gave error and then I realized why. sum(Array) turns out to be an array of (Float64,1) while the left hand side is an element of an array. Therefore, I had to do another sum on that one element array to give me a scalar. I just could not figure out the error for almost an hour and then when I did (sum(foo)) it would give me a number but S(i,j)=sum(foo) did not and gave error (I do not recall the error but I think it was something like "conversion for this type not possible"!) Of course being fresh to Julia I did not figure this out and then I realized that this occurs usually in scipy too. But you are correct, I will change the sum(sum) provided I am able to figure out how to convert the array to a number.
Also I tried plain summation without using sum just like your last section but it is was twice slower. I have that code too and I will recheck again to see how it differs from your last section. John, I will implement some of Stefan's suggestion and run the timing that you suggested. But what I am doing is very simple. I copy paste each section of the code (both Matlab and Julia) and run them and time them. I find that only that last piece which I wrote takes full 40 odd seconds. But I should find out which line takes up that time. Thanks On Wed, Jan 29, 2014 at 12:59 PM, Stefan Karpinski <[email protected]>wrote: > This sum(sum(foo,2)) business is really wasteful. Just do sum(foo) to take > the sum of foo. It's also better to extract the dimensions into individual > variables. Something like this: > > function runave1(S,A,f) > s1, s2 = size(A) > p1 = f+1 > for n = f+1:s2-f-1, m = f+1:s1-f > S[m,n+1] = S[m,n] + sum(A[m-f:m+f,n+p1]) - sum(A[m-f:m+f,n-f]) > end > S > end > > > I suspect that since Matlab forces this sum(sum(X,2)) idiom on you, it > probably detects it and automatically does the efficient thing. It's > unclear to me why you need two sum operations when the slices you're taking > are just single columns, but maybe I'm missing something here. > > Currently, taking array slices in Julia makes a copy, which is > unfortunate, but in the future they will be views. In the meantime, you > might get better performance by explicitly using views: > > function runave2(S,A,f) > s1, s2 = size(A) > p1 = f+1 > for n = f+1:s2-f-1, m = f+1:s1-f > S[m,n+1] = S[m,n] + sum(sub(A,m-f:m+f,n+p1)) - sum(sub(A,m-f:m+f,n-f)) > end > S > end > > > And, of course, there's always the nuclear option for really performance > critical code, which is to write out the summation manually: > > function runave3(S,A,f) > s1, s2 = size(A) > p1 = f+1 > for n = f+1:s2-f-1, m = f+1:s1-f > t = S[m,n] > for k = m-f:m+f; t += A[k,n+p1] - A[k,n-f]; end > S[m,n+1] = t > end > S > end > > > Not so elegant, but probably the fastest possible version. Ideally, once > array slices are views, the simpler version of the code will be essentially > equivalent to this. It will take some compiler cleverness, but it's > certainly doable. It would be interesting to hear how each of these > versions performs on your data. > > On Wed, Jan 29, 2014 at 12:30 PM, John Myles White < > [email protected]> wrote: > >> Can you show the call to @time / @elapsed so we know exactly what's being >> timed? >> >> -- John >> >> On Jan 29, 2014, at 9:28 AM, Rajn <[email protected]> wrote: >> >> > Now it takes even longer i.e., ~1 minute >> > >> > Does this make sense. Also I am running this loop only once. I do not >> understand why writing in the function form would help. I read the manual >> but they suggest writing function form for something which is used many >> times. >> > I=runave(S,A,f) >> > showim(I); >> > >> > function runave(S,A,f) >> > imsz=size(A); >> > p1=f+1; >> > for n=(f+1):(imsz[2]-f-1) >> > for m=(f+1):(imsz[1]-f) >> > >> S[m,n+1]=S[m,n]+sum(sum(A[m-f:m+f,n+p1],2))-sum(sum(A[m-f:m+f,n-f],2)); >> > end >> > end >> > S; >> > end >> > >> > Do I have to declare function parameters to speed it up. >> >> >
