onalark
So it's scale, politics, and history/momentum at this point.
[1] https://en.wikipedia.org/wiki/U.S._Army_Corps_of_Engineers_c...
Also, NumPy is implemented completely in C and Python, and makes extensive use of CPython extension hooks and knowledge of the CPython reference counting implementation, which is part of the reason why it is so hard to port to other implementations of Python.
One opportunity I do find in this article is real-world performance tests and determining the importance of constants in performance tuning. Prakop et al. did some really interesting work in 1999 on cache-oblivious algorithms, with applications to kernels such as the FFT and matrix multiplication. The work is theoretically extremely interesting, but in practice, hand-written or machine-generated algorithms continue to dominate. The most famous example of this is probably Kazushige Goto, whose hand-optimized GotoBLAS (now under open source license as OpenBLAS) still provides some of the fastest linear algebra kernels in the world.
If you're interested in learning more about how the differences between the two approaches shake out in linear algebra, I recommend "Is Cache-Oblivious DGEMM Viable?" by Gunnels et al.
A few reactions from reading through the release:
Scientist appears to be largely limited to situations where the code has no "side effects". I think this is a pretty big caveat, and it would have been helpful in the introduction/summary to see this mentioned. Similarly, I think it would be nice to point out that Scientist is a Ruby-only framework :)
You don't mention "regression test" at any point in the article, which is the language I'm most familiar with for referring to this sort of testing. How does a scientist "experiment" compare to a regression test over that block of code?
Anyway, thanks again for writing this up, I'll be thinking more about the Experiment testing pattern for my own projects.
This looks a lot like the Jupyter/IPython Notebook, which is a free and open source "scientist's notebook". If you're interested in mixing LaTeX, Markdown, and code from almost any language (Python, R, and Julia are very well-supported but there's an open kernel spec), then this might be a more appropriate tool for you to use.
The Jupyter/IPython notebook default storage format is JSON, which makes it a little more friendly for text-based version control, and also enables a static HTML view of notebooks (http://nbviewer.jupyter.org/github/ketch/teaching-numerics-w...) on GitHub.
Helen Shen wrote up a great article for Nature (http://www.nature.com/news/interactive-notebooks-sharing-the...) on how scientists are using the notebook, but it also provides a good overview of how you might use it, as well as a free interactive demo.
edit: On reviewing, I think the intent of the original blog post was to modify images in place (or at least to do it as quickly as possible with in-place filtering ok). In that case, I think my comparison is fair, since NumPy doesn't offer a faster way to do the requested operation. I didn't try out einsum, but I think Numba would outperform that as well.
Would you believe that Numba is 4 times faster for the sort of simple transformations described in the blog post?
(See aktiur's response below, some performance gains come from avoiding a copy)
As a few others have pointed out, sending a lambda function through NumPy is almost always the last thing you want to do. Unfortunately, you were in a situation where you were either going to have to do something really painful like using einsum: https://stackoverflow.com/questions/29989059/matrix-multipli... or writing your own ufunc:https://docs.scipy.org/doc/numpy-dev/user/c-info.ufunc-tutor...
I suspect that your primary limitation was the Google Compute Engine infrastructure. I'm not familiar with the limitations there, but a quick search on Google turns up a fairly limited set of libraries indeed.
I thought it would be interesting to adapt your code slightly to use Numba acceleration. Here's what it looks like:
from numba import jit
def avg_transform(image):
m, n, c = image.shape
for i in range(m):
xi = image[i]
for j in range(n):
avg = xi[j].sum()/3
xi[j][:] = avg
return image
fast_avg_transform = jit(avg_transform, nopython=True)
Re-reading your post, I suspect that einsum might actually be your cup of tea, but I really enjoy the simplicity and performance of using Numba for these sort of tasks.
Also, there are a few components on these supercomputers that you won't find on a typical workstation or cluster machine.
First, the kernel itself is frequently very small, lightweight, and much closer to what you would find on an embedded system than on a traditional desktop computer. That's because once the program on a supercomputer is loaded, the kernel's job has been mostly to get out of the way. This isn't to say that you don't find the Linux kernel on these supercomputers, it just isn't as common as you would think from reading this piece.
Fear not, many pieces of these operating systems are still open source. Here's IBM's fusedOS prototype:https://github.com/ibm-research/fusedos.
Second, many of these computers primarily only run code in C, C++, Fortran, and Python. These tend to be the only major languages in play on the HPC machines, with acceleration frameworks such as OpenMP, OpenCL, and CUDA playing major roles.
Finally, everything is glued together with MPI, a high-level (at least it was in the 90s) abstraction for scientific programming that maps down to very high-performance networks designed to help scientific codes "scale", that is, run effectively when millions of cores are simultaneously engaged.
These are beautiful machines producing important science, and the GNU/Linux operating system plays an incredibly important role in both their implementation and culture.
I'm trying to make the point that C++ occupies a weird place in the abstraction-vs-performance tradeoff.
For high abstraction and friendliness, C++ doesn't do nearly as nice of a job as MATLAB or Python in exposing high-level scientific computing operations.
For performance, C++ does have some really nice high-performance libraries (Eigen comes to mind), but for the majority of use cases, you either need to write inline assembly to get the best possible speed, or you're already reusing a fast numerical kernel from another library, and you may be better off in a higher-level language.
This is a very subjective discussion, and there are many interesting, fast, C++ libraries that do important work for science. But your claim that C++ is the best language for such work would be highly contested by many high performance computing specialists. As it is, I believe there's a pretty rough split between C, C++, and Fortran supporters, with Python continuing to gain traction.
Here's another counter-argument. MPI is used in 99% of all scientific high performance computing codes. The latest version of MPI, MPI-3, drops explicit support for C++ bindings. If C++ was so dominant, why would explicit bindings be dropped?
Many of the highest-performance scientific computing libraries are written in assembly or use inline assembly instructions. The vast majority of scientific programmers don't need to write at this level. I challenge you to write a native C++ matrix-matrix-multiply algorithm that generally outperforms the matrix-matrix-multiply available to MATLAB, Python/Cython, and Julia through their interface to the assembly-optimized OpenBLAS library: http://julialang.org/
Note that C, Fortran, and Julia are all using OpenBLAS in the rand_mat_mul benchmark, and OpenBLAS is written in C (with assembly).
The Parallella brings distributed-memory programming in, which is a very important development.
You and I strongly disagree on the meaningfulness or definition of the term supercomputer. Here's an easy definition: A supercomputer is any single, unified, computer system that is currently one of the fastest 500 in the world.
Yes, these are awesome, full-featured models, but the differences between this and a supercomputer, which costs tens of millions of dollars, requires high-density power and cooling, features multi-dimensional, low-diameter networks, and contains hundreds of thousands to millions of compute cores is... quite vast.
Let me turn this around on you. Given the availability of University preprint servers, what's the value of a service like arxiv?