Using D and std.ndslice as a Numpy Replacement

> [?] QR: Afraid I'm not sure what you're referring to.

> [+] Cholesky Decomposition: Trivial to implement

How is it that you can judge Cholesky decomposition to be trivial to implement, while being unfamiliar with QR decomposition?

Cholesky decomposition is not trivial to implement well.

> Your post comes across to me as being rather cynical - but why not be supportive of the good work that's being done?

I had a similar experience to what's happening here myself. I worked on a project that had some very specific requirements for a high performance simulation (described by a differential equation). I worked on a very specialized piece of code for a long time to optimize this differential equation solver, eventually writing a JIT compiler to generate specialized code for each simulation.

I had a D advocate come along and rather vigorously argue that D can solve this problem, without really understanding it, and missing some key requirements.

This sort of comes off the same way, "a numpy replacement". numpy is big and does a lot of stuff. Like Cholesky decomposition, many of these things have trivial implementations, and then they have good implementations. The trivial implementations will take a few hours to write. The good implementations will take a few months, or even longer. Sometimes "not good" means slow, sometimes "not good" means numerical stability issues, etc. If you don't have the right background, these things may not even be obvious.

This is a bit of a pattern with D folks. Several times (this, my previous experience, and others), I've gotten the impression that a D fan comes along, invests a bit of time (a few weeks, months) in a particular domain, and starts making comparisons with other tools that they only scratch the surface of. It can rub people the wrong way.

edit: The symbolic expression reference in this post is another example. Computer algebra is an incredibly difficult thing to implement. I don't know of anything that is harder to implement well than a computer algebra system. Yet, you link a single file under 1000 lines as partially checking the box in your list for symbolic expression manipulation. D might be a highly productive, but it's not that productive.

I occasionally rewrite Python+NumPy signal processing code in C++ for purposes of packaging and integration with native apps, so I read these examples with an eye to how they compare with typical C++, rather than with NumPy. They compare very well, and it would never have occurred to me to look into D as a possibility for this sort of code.

I'm guessing the GC might rule it out for many cases where you do signal processing in C++, but I may as well ask: what's the deployment side of things like? Can I easily build a shared library and use it from a C++ application?

> This is a bit of a pattern with D folks

Nothing D specific here. That's the pattern for every language, be it Haskell, Go, Common Lisp, JavaScript or Java ...

Numba certainly does not obviate all of these issues. I think user `tadlan` is referring to some of the compiler optimizations, like loop unrolling or fusion (e.g. noticing that two subsequent loops can be 'fused' into the subordinate execution block of one single loop). These things can offer speed-ups even beyond NumPy, and they can work even when the Python code you start out with already uses NumPy.

The thing is, which `tadlan` seems unaware of, a lot of this stuff just fails in production environments and hits corner cases that the Numba compiler does not handle (I'm talking about the first part of the Numba compiler pipeline, where it converts to Numba IR, and not yet to LLVM IR, and does things like examine the CPython bytecode to alter the representation from a CPython stack-based representation to a register-based representation that will be compatible with LLVM and ultimately with the actual machine itself). In that compilation step, the only things that are able to be handled are things that the Numba team (I used to be a member of it) explicitly support. They don't support a full-blown compiler for the entire Python language, nor even for every type of NumPy operation. That's not a knock against Numba at all -- it's a specializing compiler and obviously they need to prioritize what to support, and make longer term goals about supporting more general things. But the point still remains that you cannot just assume that if you call `jit` on any arbitrary Python code, it will always become faster. In some cases, it can even become slower.

I suspect `tadlan` is very interested in Numba and enthusiastic about knowing the taxonomy of Numba details, but it does not seem like that user has had real world experience trying to get Numba to work in production, and seeing all of the numerous buggy and missing features. I don't want to diminish anyone's enthusiasm for Numba, so it's probably best just not to engage with `tadlan` about it. That user's mind seems made up already.

"how do I know that the calculations aren't actually optimized away in all that code I can't understand? E.g. what happens if you make f0() return the means array instead of being a void function?"

If you can't understand the code, how did you know it was a void function. Please stop with the hyperbole, it's not adding anything to the discussion.

Updated code:

    import std.range : iota;
    import std.array : array;
    import std.algorithm;
    import std.experimental.ndslice;
    import std.datetime;
    import std.conv : to;
    import std.stdio;

    enum testCount = 10_000;

    auto f0() {
        auto means = 10_000_000.iota
            .sliced(1000, 10000)
            .transposed
            .map!(r => sum(r) / r.length)
            .array;
        return means;
    }

    void main() {
        auto r = benchmark!(f0)(testCount);
        auto f0Result = to!Duration(r[0] / testCount);
        f0Result.writeln;
    }

Results:

    Python: 14.1 msec
    D:      41   μs
    D is 343.9x faster

"Given the same number of lines (or characters), are you sure you can't write equally fast Python code?"

IMO program size is an almost meaningless statistic outside of code golf challenges. LOC is not an indicative measure of code readability, usefulness, or organization.

For example, your Fortran code was 11 lines while the D function (with the return) is eight lines.

"I tested the Fortran version on a slightly slower computer (got 14.7 msec on the Python version). Fortran runs at 0.7 μs, i.e. > 55x faster than D, with less code that's more readable to boot."

Just goes to show why Fortran is still used in a lot of scientific areas.

But two things:

1. Fortran is a much simpler language than D or Python and is much harder to do multipurpose work in it (so I'm told from Fortran programmers, I don't know Fortran myself). So when your program needs to do anything other than number crunching, it's normally done in a separate language. Using D you can have everything in one code base.

2. This article was about Numpy and std.ndslice because those are two areas that I know about and Numpy is a very popular library. Bringing up Fortran's speed here is like commenting on how much faster C++ is in a thread about Ruby.

Also, readability is a subjective idea; I believe the D code is more readable than the Fortran you wrote. Different strokes.