Cram embed_ref=2 and huffman optimisation

[jkbonfield]

• Put an upper limit on the number of cram huffman symbols to prevent malformed data consuming a lot of ram and cpu.
• Prevent excessive memory and/or cpu usage in embed_ref=2 for extremely sparse data or huge deletion / ref-skip CIGAR ops.
• Plus a trivial compiler warning fix for fuzz builds.

[daviesrob]

This is over-complicated, especially as hist_t is now uint8_t.

Suggested change

	#define MAX_HIST_T (uint32_t)((1LL << (8*sizeof(hist_t)))-1)
	#define MAX_HIST_T UINT8_MAX

[jkbonfield]

Agreed. It's better than the alternative I was thinking of anyway with const uint32_t instead.

[daviesrob]

Could we remove the +10 on both sides of the comparison?

[jkbonfield]

As it stood before, it passed the test on tiny containers of 3 records and switched to multi-seq mode.

It's lost in the mists of time[1] as to why I had the +10 there before, but it was wrong anyway. I think I was trying to do a small sample correction, but that's generally (a+const)/(b+const) rather than one-sided. The small constant basically means we don't take these decisions to switch mode unless there is a significant number of alignment records, which I wish to keep intact.

[1] Ie it's maybe in io_lib's history somewhere, but it predates the code arriving in htslib.

[jkbonfield]

Found it. It arrived in 2013. :-)

jkbonfield/io_lib@ecec3e6

There was a lot of code though as it was the creation of the multi_seq mode so the finer nuances weren't documented,

[daviesrob]

It may have made sense when the +10 was only on one side, but now it's on both and c->curr_rec+10 < c->max_rec/4+10 is the same as c->curr_rec < c->max_rec/4 (ignoring the possibility of overflow).

[jkbonfield]

Doh sorry, implementation failure. It's meant to be (c->curr_rec+10) < (c->max_rec+10)/4 I think for a small sample size correction. What I had before vaguely made sense if it was on the left side (or a -10 on the right obviously), but I'm not sure quite what I was thinking as it was 13 years ago.

Or put another way (if it was floating point): (c->curr_rec+10.0)/(c->max_rec+10.0) < 0.25.

[jkbonfield]

Basically compare this:

$ perl -e '$N=20;for ($i=0;$i<$N;$i+=$N/10) {print $i/$N,"\t",($i+10)/($N
+10),"\n"}'
0	0.333333333333333
0.1	0.4
0.2	0.466666666666667
0.3	0.533333333333333
0.4	0.6
0.5	0.666666666666667
0.6	0.733333333333333
0.7	0.8
0.8	0.866666666666667
0.9	0.933333333333333

$ perl -e '$N=2000;for ($i=0;$i<$N;$i+=$N/10) {print $i/$N,"\t",($i+10)/($N+10),"\n"}'
0	0.00497512437810945
0.1	0.104477611940299
0.2	0.203980099502488
0.3	0.303482587064677
0.4	0.402985074626866
0.5	0.502487562189055
0.6	0.601990049751244
0.7	0.701492537313433
0.8	0.800995024875622
0.9	0.900497512437811

If we wanted to trigger behaviour when we drop below, say, 50% capacity then on a large block that's basically curr/max at 0.5. On a small block we need a more extreme ratio of 0.2ish, to compensate for error from a low sample size.

The +10 is entirely arbitrary of course, but the principle is sound. Unlike the original (or amended) code!

[daviesrob]

Well, that makes more sense, although the rewritten version is really just c->curr_rec+7.5 < c->max_rec/4. And to avoid the floating point, you'd probably go for either 7 or 8 on the left.

[daviesrob]

Apologies for incorrect original version of that, due to mathematical ineptitude on my part...

[jkbonfield]

I now have c->curr_rec+8 < (c->max_rec+8)/4 which can only succeed when max_rec >= 32. I'm not sure if this is the correct ratio, but it's sufficient I think to act as a starting point.

All of this is heuristics anyway in how best to optimise the CRAM file and when to switch into multi-ref mode.
The intention was simply if you're consistently producing CRAM containers which are a 1/4tr of their requested size, in terms of the number of alignments you've put into them containing the same reference, then switch out of single-ref mode into multi-ref mode (and use the RI data series to indicate which reference is which instead of a constant value in the slice header).

The purpose was for unsorted data basically, where we were originally creating thousands of small containers with 1 or 2 reads in as the chances of have many reads for the same reference together is slim. We just don't want to switch to this mode inappropriately, but the bar is deliberately quite high.

[jkbonfield]

Agreed on it basically being able to be simplified, but it obscures the derivation. (X+c)/(N+c) is a standard mechanism for small sample size correction in significance tests, but factoring it out just hides this.

It'll all be reordered by the compiler anyway I expect so let's stick with something clearer. Feel free to add a comment explaining the origins of it if you wish, although it's documented here for anyone who wishes to delve into git blame.