Fix octopus binary "invalid instruction" (-march too new)

[holtgrewe]

• I have read the guidelines for bioconda recipes.
• This PR adds a new recipe.
• AFAIK, this recipe is directly relevant to the biological sciences (otherwise, please submit to the more general purpose conda-forge channel).
• This PR updates an existing recipe.
• This PR does something else (explain below).

@chapmanb @kyleabeauchamp What do you think?

[chapmanb]

Manuel;
Thanks for looking at this. Would you be able to explain more about the motivation and consequences of this? Not being a C++ expert, it seems like this would be less portable by requiring intel sandy bridge chips on any machines downloading and using these binaries. I don't know what the defaults do without this, so could use a little handholding about the approach you're taking. Thanks again.

[holtgrewe]

@chapmanb Ah, sorry for not giving much context here. We are hitting an issue here with that the binaries for octopus actually require a Haswell architecture (one or two generations after Sandybridge) and we get an "invalid instruction" error on our Sandybridge cluster nodes. At another point, I asked about policies of this but there really was none somehow.

I think #2354 contains some relevant discussion and links. Maybe @kyleabeauchamp can give some insights into this. I would be fine with "going back" further or even somehow shipping multiple architectures at the same time.

Right now, the package is broken for us and I could imagine that this is the case for people who have even older cluster nodes (e.g., high memory nodes tend to be kept alive longer than the average compute node as they are pricy).

[chapmanb]

Manuel -- thanks much for the context. I'd definitely be keen to make this as portable as possible. I don't know enough about this to comment usefully, but it would be helpful to have @dancooke in on the conversation. Dan, we're working on making the octopus binary more portable. Do you know the potential impact of compiling against older or multiple architectures on runtime? Any tips/tricks for us doing this in the best way?

[epruesse]

Recipe's should not set -march or -mtune at all. They should actually take care to remove any such flags from upstream build systems. (I'm surprised the CB3 compilers don't flat out override those flags.)

The main argument here is that if someone needed to eek out those extra 3%, they could take the parts of Bioconda they need, adjust the conda_build_config.yaml, and rebuild into their own channel. If the recipe messes with those flags, that might not work as desired.

The someone might be us at some point. Or just you. Or perhaps a new project called Biogentoo. In any case, the default for Bioconda at this point needs to be "works everywhere".

just my $0.02

[holtgrewe]

@epruesse I agree. But how do we achieve this in Bioconda for this package in a short timeframe?

[epruesse]

@holtgrewe Just add a patch removing the -march=native bit from src/CMakeLists.txt and open an issue upstream, notifying them that their Release builds code that won't run on a CPU older than the build machine.

[dancooke]

@chapmanb Thanks for looking into this. I can't say that I've done much benchmarking across different architectures. However, I did find that compiling on one of our centres cluster nodes that supports AVX-512 did speed up runtimes considerably - most noticeably for somatic calling - compared with nodes that only support AVX or AVX2. I could imagine this is due to better vectorisation in the maths heavy methods used for somatic calling.

[epruesse]

@dancooke All 64 bit CPUs have SSE2, so we always get that. Depending on what's being done, SSE4.2 brings significant extra performance, while AVX and AVX2 aren't worth it in my experience. AVX-512 again brings a lot, but it's actually a collection of extensions of which CPUs only implement some.

Options:

• follow the RAxML approach and compile a bunch of variants of the binary.
• go above that and have a wrapper script selecting the right one.
• go the vsearch route and do this inside of your binary with hand optimized SIMD instructions (although vsearch doesn't get this 100% right and breaks on a small subset of the Travis nodes - testing is really hard).
• build the compute kernels in a separate DSO, build that multiple times for a range of -march targets and dlopen the right one at run time
• create a channel, or multiple channels, of your own, in which you place copies of the packages built for specific architectures. You can then layer the appropriate channel above Bioconda to override packages with custom optimized ones.
  Bear in mind though, that this needs a separate miniconda folder for each CPU architecture.

[chapmanb]

Thanks for all this helpful discussion. My take here is that we should merge Manuel's fix. The upstream approach provides some value on specific architectures if pre-compiling so will likely stay in place. I don't see anyone jumping to try and maintain multiple optimized versions in bioconda right now, so the fix here will give us better compatibility without an obvious runtime issue. What does everyone else think?

[holtgrewe]

I think we need @epruesse's suggestion of patching the CMakeLists.txt. I can look into this tomorrow.

[epruesse]

There you go. Let's see if that works.

[lh3]

Intel developers suggest that the right strategy is to use the so-called CPU dispatch. The idea is to compile different versions of the same functionality at the compile time. During run time, the program dynamically chooses the fastest version based on the host CPU. GATK uses this approach.

ksw2_dispatch.c from minimap2 shows a working example. It prefers sse4.1 if present, or falls back to sse2 otherwise.

[epruesse]

@lh3 Thanks for the links!

In the end, dynamically selecting code paths suited to the CPUs capabilities at run time is the way to go.

It can happen at whole-binary level (e.g. raxml-AVX - though it's missing a dispatcher), or library level (e.g. Intel MKL, where libmkl_core contains a dispatcher and libmkl_avx2 and friends are optimized), or function level (e.g. minimap2 or vsearch).

The higher levels are easier if you want to let the compiler optimize for features, the lower if you want to do it by hand.

The catch with hand optimized SIMD code, such as vsearch and minmap2 are using, is that it's not just hard to write, it's also hard to test without having a large selection of computers with different CPU generations on hand. (We have an issue here regarding a vsearch problem on some Travis instances...)

The catch with letting the compiler do the optimization is that the tree vectorizers are somewhat fickle with regard to the way your loops are written. So you have to monitor with diagnostic flags that your code is actually still vectorized. And, as happened to DADA2, you may find that users with older compilers suddenly loose a lot of speed (Bioconda builds with GCC 4.8 were orders of magnitude slower).

No free lunch anymore... :)

[lh3]

When vectorized, some algorithms take a distinct form. Automatically vectorizing such algorithms is beyond the capability of compilers. Hand written intrinsics is the only option. Debugging CPU dispatch is actually easy. Because dispatching is dynamic, you can implement a command-line option to target a particular instruction set at run time. You only need one recent computer to test all older architectures.

CPU dispatch does pose a challenge to conda because it may require fairly recent CPU and compilers. Purely from an enduser point of view, the best solution is to provide a portable binary precompiled on a recent machine (GATK follows this route). The conda recipe only downloads the binary, not compiling from the source. This way, all users get the maximal performance on their machines. From a developer point of view, however, this can be very complicated to implement.

[dancooke]

@lh3 This is an interesting approach, but I would guess that the number of routines where vectorisation results in a meaningful difference to overall runtime, and that the compiler can't automatically vectorise, is small for most programs. In Octopus, the only algorithm that uses hand-coded SIMD instructions is the pair hMM, which uses SSE2. In this routine, the size of the SIMD registers determines the band-size of the DP table, so using a different SSE instruction set actually changes the behaviour of the algorithm (by allowing bigger gaps). Of course, the SIMD code could be modified to keep functionality equivalent to the largest register variant, but then this penalises the smaller register variants by enforcing a larger band-size. The band-size should be defined by the program or user, not the system architecture, so we'd also need to write routines for each possible band-size and select the correct implementation based on that and the system architecture.

On the other hand, compilers are very good at vectorising other methods, such as inner products, that account for much of the computation time in Octopus (e.g. in expectation-maximisation). Writing custom dispatch-based SIMD code for these methods would likely result in little/no performance improvement, but a maintenance headache.

So ultimately, I think that the high-level approaches suggested by @epruesse are the way to go. The multiple-binary 'wrapper script' approach seems a good option. @epruesse, I presume this means the script selects the correct binary on installation, rather than at runtime?

Thanks everyone for all their help and suggestions on this issue.

[lh3]

On paired-HMM, you are explicitly using SSE2. Keeping it that way requires no changes to the code. On inner products etc, you can create different versions of sdot() and choose the best version at runtime. Another option is to use a BLAS library and let the library deal with the rest. That is one extra dependency, though. Note that the point of CPU dispatch is to let users conveniently and consistently get the best performance out of our tools. CPU dispatch does make development more difficult. That is the price we developers pay for better user experiences.

If you go for the wrapper approach, the wrapper should choose the right version at runtime. In a cluster environment, we install a tool on one machine and run it on other machines that may have different specs. The launcher script has to choose the right version at runtime, or the tool may cause the invalid instruction error (if the install machine has a higher spec) or underperform (if the install machine has a lower spec).

By the way, back to this PR, I fully agree that we should remove "-march=native". CPU dispatch is more for future considerations.

[epruesse]

I presume this means the script selects the correct binary on installation, rather than at runtime?

@dancooke You can't assume that the CPU features available at installation time are also available at run time. Small clusters in which new generations of compute nodes are added over time are probably pretty common. So I'd definitely detect the CPU at run time: the cost in storage and CPU cycles is negligible, whereas explaining to each new Bioconda user how CPUs are different is anything but.

[epruesse]

(Build failed because conda.anaconda.org had a 5xx error - let's retry that)

[epruesse]

[off topic]

When vectorized, some algorithms take a distinct form. Automatically vectorizing such algorithms is beyond the capability of compilers. Hand written intrinsics is the only option.

@lh3 My experience is very limited, so not contesting your claim here. When I was working on ARB's parsimony implementation, I found it possible to write that "distinct form" in C, allowing the compiler to vectorize while avoiding intrinsics. The code ended up looking rather unusual, and required checks in the build system to prevent minor edits from breaking vectorization. The hope was that this would be more future proof than assembly style direct use of intrinsics.

Debugging CPU dispatch is actually easy. Because dispatching is dynamic, you can implement a command-line option to target a particular instruction set at run time. You only need one recent computer to test all older architectures.

While that allows you test each implementation, it doesn't allow you to test whether an implementation works on a CPU with a smaller instruction set. You could still have an AVX2 instruction in your Pentium4 implementation and wouldn't catch it (in reality it's going to be more subtle than that of course). I think that's what's causing the vsearch issue.

I wish there was some kind of test framework that allowed checking this. Never found anything though.

CPU dispatch does pose a challenge to conda because it may require fairly recent CPU and compilers. Purely from an enduser point of view, the best solution is to provide a portable binary precompiled on a recent machine (GATK follows this route). The conda recipe only downloads the binary, not compiling from the source. This way, all users get the maximal performance on their machines. From a developer point of view, however, this can be very complicated to implement.

Why does the compile host matter? You may not be able to run tests, but "cross compiling" for more capable CPUs should only be a function of the compiler.

BTW: There is a WG21 proposal for extending C++ with standard templates abstracting from SIMD intrinsics (P0214R9, partial implementations in dimsum, vc, pik/simd).

[lh3]

We are talking about slightly different things, but anyway I agree with you. If you know what assembly or intrinsics-based code should look like, you can write unusual but pure C code and let the compiler vectorize for you. I know some Intel developers are doing this. However, their implementation at the time only worked with Intel compiler.

On debugging, you can disable advanced instruction sets at the compile time, such that calling those instructions leads to a compiling error.

Why does the compile host matter? You may not be able to run tests, but "cross compiling" for more capable CPUs should only be a function of the compiler.

That is host dependent. Gcc on several of my machines are recent enough but they don't have advanced instruction sets built in.

[epruesse]

@bgruening Can you take care of the "broken" build?