Is hts_os.h part of the public API?

[jmarshall]

Is it really intended that the recently added htslib/hts_os.h is part of HTSlib's public API?

If it is an HTSlib-internal thing, would it not be better as os_internal.h or so, in the top-level directory and in particular outside the public header directory.

[jkbonfield]

Given that samtools depends on it, it has to be public. The intention is to provide missing OS functionality so that htslib can rationalise code that expects a POSIX interface. It's not ideal, but the alternative is to duplicate such code in both packages, which also isn't ideal.

[jmarshall]

Good luck with that. It needs more #ifdef trickery to avoid screwing up platforms that do actually get this right.

[jkbonfield]

Care to comment what then? It would have been good to know about this in the month the PR was waiting to be merged too!

[jmarshall]

I did a quick search for “don't declare standard library functions yourself, use the system headers” but didn't immediately find a good reference in e.g. the C FAQ. In any case, (1) C99 §7.1.3/2; (2) this is a long-standing good practice that I know you're more than familiar with; (3) you just had a refresher of this in cram/os.h.

So for example those declarations of [sdl]rand48() ought to be omitted on platforms where <stdlib.h> has already provided them, possibly doing so as macros or with interesting annotations — in which case hts_os.h's versions would conflict with the system ones. This could be via configure but that's extremely fragile for a public header, or hopefully something like #ifdef _POSIX_VERSION or similar would do the trick.

Care to name the relevant samtools PR?

Presumably you have colleagues who are still in a position to review PRs as part of their paid jobs 😄

[jmarshall]

This replacing system functions malarky comes down to the same considerations that previously resulted in HTSlib introducing its own hts_md5_init() et al MD5-equivalent functions, rather than providing MD5_Init() et al with their usual names and attempting the gymnastics required to avoid colliding with system-provided functions (here “merely” declarations, mostly) with the same names.

[jkbonfield]

Agreed the extern definitions of drand et al ought to be ifdef guarded and checked via autoconf for consistency and the fact they're not is simply jjust me forgetting as the others have various ifdefs around. As you saw, I added HAVE_DRAND48 already to the Makefile. However they are the POSIX prototypes so haven't yet caused issues on any POSIX compliant systems.

It's possible they also break building under _WIN32 with other compilers as the existing guards aren't checking HAVE_SRANDOM etc, but we know MSVC has many more issues and it doesn't work there, so we'll cross that bridge when we come to it. Ideally it'd all be done by autoconf, but I don't wish to go down that road while we have a schizophrenic build system that is both autoconf and not autoconf at the same time as it simply doubles the work load for handling this sort of thing.

As for the review - it was checked by more than one person, yourself included. The related samtools PR is samtools/samtools#678. You'll see there are quite a few samtools things that use the rand48 functions, such as bam_sort.c. Similarly in bcftools too. We could change all of those to use another random number generator instead or emulate 48 bit precision by combining multiple 16-bit or 32-bit random number generators, but frankly we're targetting a basic POSIX system here so just including a sufficient compatibility layer to get it working is, hopefully, sufficient. Time will tell whether this causes issues with other packages, but it's meant to be limited to Windows systems and we can address these as and when they arise.

[jmarshall]

Indeed, the rand48 functions are used extensively in samtools but (rather to my surprise) currently not at all in bcftools. In htslib, they are used only in ks_shuffle(), which really needs only an integer RNG with fewer bits. So my plan was to change the relevant ks_shuffle() line to j = KSORT_RAND(i) where KSORT_RAND is a macro defined in terms of drand48() if obviously available otherwise rand(). Then any rand48 heroics would be limited to samtools.

As you say, HTSlib primarily targets POSIX. The level of Windows support is a policy decision, and sensible levels could be:

1. doesn't even build out of the box (which is what's historically been the case)
2. builds out of the box, and mostly works (possibly not identically to other platforms)
3. builds out of the box, passes the entire test suite, works just like on POSIX

I was aiming at (2), in which case a simple stubbing out of the rand48 functions for — say — rand() equivalents on Windows would suffice. So a sufficient compatibility layer would be a set of one-liners for the rand48 functions and it wouldn't matter how many of samtools/htslib/bcftools it was duplicated across.

Providing an implementation of the real rand48 functions for Windows is a lot more than the minimal sufficient solution. Is it that the samtools maintainers have made a conscious decision to aim for (3)?

As for the review - it was checked by more than one person, yourself included.

For the record, I comment on things that catch my eye or that affect parts of HTSlib that I am interested in, and hopefully that is useful. That does not constitute a full review, and as it is no longer my day job I don't have the time to provide full reviews.

[jkbonfield]

IMO, and I think of the team, is that method 2 isn't helpful as you don't really know whether or not it works if the test harness doesn't pass, and you can't spot regressions (of which there appear to have been many). So it comes down to don't support at all vs support properly.

Given we now have AppVeyor for CI testing, method 3 seems the best, plus we now have in house testing & development on Windows too. We do have users who run the code on Windows too (some instrument manufacturers like to do this as their instrumentation software is already on windows, so a native binary avoids complex emulation / VM infrastructures).

Clearly it's all a work in progress and I have no doubt we'll find subtleties here and there, but that's the nature of the develop branch. I'll tidy up the ifdefs though.

[jkbonfield]

Thinking more on this, the other issue is potential speed. Using rand()%val looks to be around 30% slower than drand48()*val due to the slow division, in addition to needing more work to cope with larger ranges of val. I've no idea how significant that is though. I'd be very surprised if it's the bottleneck in existing code that uses ks_shuffle (eg variant calling).

I made a call on implementation details and the absence of any timely feedback meant it was accepted, but we should (and will) discuss the specifics of our non-POSIX strategy within our weekly meetings though rather than simply on the PR where it gets lost in a sea of other tweaks.

[lh3]

On pseudorandom number generator (PRNG), another possibility is to implement one by ourselves. This would improve randomness and guarantee identical random sequence on different platforms. There are a few options on how to implement a PRNG.

1. With LCG. This is used by glibc rand() and lrand48(). Its randomness is not that good, though.

2. With splitmix64. It is slightly slower than LCG according to the Big Crush paper [PDF], but produces better random sequences and is fairly simple (see appendix). I think this is the default PRNG in Java.

3. With more advanced PRNGs such as the xoroshiro+ family or the PCG family.

My preference is splitmix64. Option 3 is also good to me. As for how to convert uint64_t to double, the xoroshiro+ website gives a fast solution (see below).

Appendix:

uint64_t splitmix64(uint64_t x)
{
  uint64_t z = (x += 0x9E3779B97F4A7C15ULL);
  z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
  z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
  return z ^ (z >> 31);
}
double uint64_to_double(uint64_t x)
{
  const union { uint64_t i; double d; } u = { .i = UINT64_C(0x3FF) << 52 | x >> 12 };
  return u.d - 1.0;
}

[jmarshall]

If the samtools maintainers have decided to have HTSlib provide something of a POSIX-compatibility layer, very laudable goal I'm sure and good luck with that. To avoid future collisions and headaches I would suggest that HTSlib's binary libraries not provide the exact same symbols that system libraries or other libraries are or ought to be providing. So for these *rand48 functions you might have libhts.so provide hts_drand48 et al symbols (on some platforms) and enable standard drand48()-etc using code on those platforms with

// hts_os.h
#include <unistd.h>
#ifndef _POSIX_VERSION  // or similar to restrict to deficient platforms
extern double hts_drand48(void);
static inline double drand48(void) { return hts_drand48(); }
// etc (or possibly even better via #defines rather than inline functions)
#endif

I would reiterate a recommendation to not have things that are installed as source (i.e., public headers) depend on the results of configure. Needing to point your compiler at different headers for different configurations is painful, as e.g. libcurl has finally managed to stop doing.

[jkbonfield]

I would reiterate a recommendation to not have things that are installed as source (i.e., public headers) depend on the results of configure.

That's actually the original reason I omitted putting ifdef HAVE_DRAND around the declarations in hts_os.h. Rob did suggest though having hts_os.h.in and generating it locally so all autoconf shenanigans are done before compile (and hence before install). Infact this is what io_lib does with its os.h.

[daviesrob]

Rob did suggest though having hts_os.h.in and generating it locally so all autoconf shenanigans are done before compile (and hence before install). Infact this is what io_lib does with its os.h.

I'm going off that idea now, having read the curl posting. We should have hts_drand48() etc., though.

I'll look to see how some other projects solve this.

[jmarshall]

That's precisely what I'm strongly recommending not doing.

Please do not make it so that the contents of any public headers change due to configure tests. Doing so makes life harder for packagers, cross-compilers, and other users.

Alternatives for this *rand48() stuff is to wrap it in #if defined _WIN32 (tweaked so that e.g. Cygwin escapes the net if necessary) if you don't mind highlighting the Windows lossage, or perhaps something based on #include <unistd.h> #ifndef _POSIX_VERSION if that doesn't catch any innocents. The point is to delay the decision until when third-party application code is being compiled, not hard-coding it sooner.

[jkbonfield]

So by and large, you're proposing more or less what we have now but with an additional guard (be it WIN32,_ POSIX_VERSION or similar)? That seems reasonable and a minimal tweak. I do want to try and sanitise systems so htslib works across the board, and it's simply less effort than renaming functions or implementing our own versions of everything just for the sake of one slightly borked setup. That feels like tail wagging the dog.

I initially thought of adding the HAVE_DRAND48 guards around it, as suggested, but then rapidly realised it just doesn't work and indeed I'd been down that road before and backed out of it sharpish. Adding those sorts of guards to a public header means all users of that header then have the obligation of checking for drand48, eg in autoconf. That's simply too onerous, and largely pointless as if you've gone to those extents then you're most of the way there to avoiding it too.

Hence my hack of just defining them always and seeing if (pretty unlikely tbh) it ever fails. Yes they could be implemented purely in macros and yes it technically is undefined to redeclare them in this manner, but in practice until it happens I didn't really deem it a problem worth solving. It hasn't yet broken anything, but I'm all too aware we still don't have a proper build farm setup for testing on all our whacky OSes, only the mainstream ones. :-(.

The _POSIX_VERSION check seems a reasonable middle ground as it's going to be sufficient to spot machines claiming to be posix, which by definition should then have their own implementation. It's basically like the HAVE_DRAND but without needing extra autoconf / makefile hackery. :)