As usual, it depends on how much you value speed vs simplicity. You can go
farther towards the simplicity direction by noting that, in the matlab-
inspired version,
broadcast(+, a, b)
is equivalent to
a .+ b
which is a much nicer syntax than is available in Matlab.
The Python version looks more tuned for speed, but the fact that you're
allocating so much memory is an indication there's much more you can do. The
manual section that Kevin pointed out will help. Basically, you want to
replace operations like
d[1:4:n, 1:4:n] = b - d[1:4:n, 1:4:n]
with something like
subtract_d_from_b!(d, 1:4, 1:4, b)
and write that function in devectorized form.
--Tim
On Monday, August 25, 2014 06:38:11 AM Phillip Berndt wrote:
> Hi julia-users,
>
> I've recently stumbled over Julia and wanted to give it a try.
>
> To assess it's speed, I've implemented another micro-benchmark, namely a
> version of Matlab's magic() function that generates magic squares. Since I
> have no experience writing optimal Julia code, I started off with literal
> translations of two different implementations - Matlab's and the one from
> magic_square.py from PyPy, which is an optimized version for NumPy. I then
> timed the calculation of all magic squares from N=3 to N=1000. The table
> from Julia's homepage suggests that in most cases, it is significantly
> faster than Python and Matlab. In my case, it's significantly slower, which
> is somehow disappointing ;) My question now is:
>
> Can the implementation be optimized to outperform the other two?
>
> *The times:*
>
> Julia, Matlab version: elapsed time: 18.495374216 seconds (13404087428
> bytes allocated, 12.54% gc time)
> Julia, Python version: elapsed time: 8.107275449 seconds (13532473792 bytes
> allocated, 26.99% gc time)
> Matlab: Elapsed time is 4.994960 seconds.
> Python: 1 loops, best of 3: 2.09 s per loop
>
> My test machine is a 4 Core i7-4600 Notebook with 2.1 GHz and 8 GiB RAM,
> running a current Linux Mint and Julia 0.3 stable. To be fair, Python does
> not seem to gc during this loop (disabling gc doesn't alter the time here),
> so one should compare with 8.1 s * (1.-.2699) = 5.91 s for Julia. That's
> still much slower than Python. (By the way, even Octave only needs 4.46
> seconds.) If I translate the matrices in magic_python to account for
> column-major storage, the execution time does not significantly improve.
>
> *The code:*
>
> Matlab: tic; arrayfun(@magic, 3:1000, 'UniformOutput', false); toc
> IPython: import magic_square; %timeit [ magic_square.magic(x) for x in
> range(3, 1001) ];
> Julia: I've uploaded the code to a Gist at
> https://gist.github.com/phillipberndt/2db94bf5e0c16161dedc and will paste a
> copy below this post.
>
>
> Cheers,
> Phillip
>
>
> function magic_matlab(n::Int64)
> # Works exactly as Matlab's magic.m
>
> if n % 2 == 1
> p = (1:n)
> M = n * mod(broadcast(+, p', p - div(n+3, 2)), n) +
> mod(broadcast(+, p', 2p - 2), n) + 1
> return M
> elseif n % 4 == 0
> J = div([1:n] % 4, 2)
> K = J' .== J
> M = broadcast(+, [1:n:(n*n)]', [0:n-1])
> M[K] = n^2 + 1 - M[K]
> return M
> else
> p = div(n, 2)
> M = magic_matlab(p)
> M = [M M+2p^2; M+3p^2 M+p^2]
> if n == 2
> return M
> end
> i = (1:p)
> k = (n-2)/4
> j = convert(Array{Int}, [(1:k); ((n-k+2):n)])
> M[[i; i+p],j] = M[[i+p; i],j]
> i = k+1
> j = [1; i]
> M[[i; i+p],j] = M[[i+p; i],j]
> return M
> end
> end
> @vectorize_1arg Int magic_matlab
>
> function magic_python(n::Int64)
> # Works exactly as magic_square.py (from pypy)
>
> if n % 2 == 1
> m = (n >> 1) + 1
> b = n^2 + 1
>
> M = reshape(repmat(1:n:b-n, 1, n+2)[m:end-m], n+1, n)[2:end, :] +
> reshape(repmat(0:(n-1), 1, n+2), n+2, n)[2:end-1, :]'
> return M
> elseif n % 4 == 0
> b = n^2 + 1
> d = reshape(1:b-1, n, n)
>
> d[1:4:n, 1:4:n] = b - d[1:4:n, 1:4:n]
> d[1:4:n, 4:4:n] = b - d[1:4:n, 4:4:n]
> d[4:4:n, 1:4:n] = b - d[4:4:n, 1:4:n]
> d[4:4:n, 4:4:n] = b - d[4:4:n, 4:4:n]
> d[2:4:n, 2:4:n] = b - d[2:4:n, 2:4:n]
> d[2:4:n, 3:4:n] = b - d[2:4:n, 3:4:n]
> d[3:4:n, 2:4:n] = b - d[3:4:n, 2:4:n]
> d[3:4:n, 3:4:n] = b - d[3:4:n, 3:4:n]
>
> return d
> else
> m = n >> 1
> k = m >> 1
> b = m^2
>
> d = repmat(magic_python(m), 2, 2)
>
> d[1:m, 1:k] += 3*b
> d[1+m:end, 1+k:m] += 3*b
> d[1+k, 1+k] += 3*b
> d[1+k, 1] -= 3*b
> d[1+m+k, 1] += 3*b
> d[1+m+k, 1+k] -= 3*b
> d[1:m,1+m:n-k+1] += b+b
> d[1+m:end, 1+m:n-k+1] += b
> d[1:m, 1+n-k+1:end] += b
> d[1+m:end, 1+n-k+1:end] += b+b
>
> return d
> end
> end
> @vectorize_1arg Int magic_python
>
> print("Matlab version: ")
> @time magic_matlab(3:1000)
>
> print("Python version: ")
> @time magic_python(3:1000)
