Oops I forgot I'd set some bcftools benchmarks off too. I meant to put those in the commit message, but no matter.

time ./bcftools.gl mpileup --fasta-ref $HREF38 /var/tmp/novaseq-10m.bam | egrep -v '^#' | md5sum
[mpileup] 1 samples in 1 input files
[mpileup] maximum number of reads per input file set to -d 250
ce10e64e4de9623e8640d82730d2d2dd  -

real	8m12.121s
user	8m29.351s
sys	0m19.585s

time ./bcftools.dev mpileup --fasta-ref $HREF38 /var/tmp/novaseq-10m.bam | egrep -v '^#' | md5sum
[mpileup] 1 samples in 1 input files
[mpileup] maximum number of reads per input file set to -d 250
ce10e64e4de9623e8640d82730d2d2dd  -

real	9m52.420s
user	10m9.672s
sys	0m19.515s

That makes it 20% faster. It's a better less of a speedup than samtools as bcftools is doing other work too. I'll take a look at that next infact. :-)

Note if we're willing to accept non-exact matching results for a further ~5% mpileup speed up, it's possible to ditch that log and replace it with an approximation that uses the IEEE fp logic to pull out the exponent. Eg:

// Tylor (deg 3) implementation of the log2(x):                                 
// https://www.flipcode.com/archives/Fast_log_Function.shtml                    
static inline double fast_log2 (double val)                                            
{                                                                               
   register int64_t *const     exp_ptr = ((int64_t*)&val);                      
   register int64_t            x = *exp_ptr;                                    
   register const int      log_2 = ((x >> 52) & 2047) - 1024;                   
   x &= ~(2047LL << 52);                                                        
   x += 1023LL << 52;                                                           
   *exp_ptr = x;                                                                
                                                                                
   val = ((-1.0f/3) * val + 2) * val - 2.0f/3;                                  
                                                                                
   return val + log_2;                                                          
}                                                                               

We could maybe use this elsewhere in htslib infact.
The differences are slight. 4353 of 36994133 qualities differed in rounding, implying 0.011% changed.

It does however make things differ. (Marginally more than fixing that .499 vs .5 double to int rounding.)

According to this SO question, the more accurate log2 is this (single-precision):

static inline float mg_log2(float x)
{
	union { float f; uint32_t i; } z = { x };
	float log_2 = ((z.i >> 23) & 255) - 128;
	z.i &= ~(255 << 23);
	z.i += 127 << 23;
	log_2 += (-0.34484843f * z.f + 2.02466578f) * z.f - 0.67487759f;
	return log_2;
}

For integers between 0–100000, it is indeed more accurate. The SO page also points to this implementation, it is even more accurate for integers in that range.

As I've got your attention Heng, I'd like your thoughts on this alternative approach.

A considerably faster approximation is to simple avoid calling BAQ except where we think it matters. I've been experimenting with putting in the BAQ call to the mpileup loop so as we add new reads to mpileup it scans the CIGAR for indels. Then on each column iteration, if any read in the current pileup contains an indel it computes BAQ (if not already done) on all reads in the pileup (potentially this can be reduced further).

The net effect is a considerable reduction to the number of BAQ calls, making bcftools mpileup run at approx double the speed. So far my assessment with SynDip at 45x and a 15x subsample is this is actually more accurate than calling BAQ everywhere, as it sometimes harms SNP calls.

It's still a work in progress, but if this works then likely the small gain by using a faster log will become somewhat redudant.

BAQ played its role in early days but I think it is outdated and its use should be minimized.

The ideal would clearly be local MSA to produce allele consensi, then BAQ becomes redundant. However I think it still has a use until such a thing arrives. Minimising it to only run on cases where we think the pair-wise alignment may not represent the full sample(s) MSA is clearly both beneficial for speed, but also accuracy as we're not reducing qualities in places where the alignment was already correct.