Parallel cram2bam

[jkbonfield]

The main thread does a lot of work in cram_to_bam and bam_set1. Io_lib's model for this was to put this code into the worker threads. I thought I'd done that here too (I have for the writer already), but apparently not. This PR fixes that.

A profile of the main thread:

$ ./test/test_view.dev -B -@32  ~/scratch/data/novaseq.10m.cram & perf record -t $! 2>/dev/null;perf report -n 2>/dev/null| head -20 | egrep -v '^#'
[1] 3634487
[1]+  Done                    ./test/test_view.dev -B -@32 ~/scratch/data/novaseq.10m.cram
    49.75%          2100  test_view.dev  test_view.dev      [.] bam_set1
     9.51%           407  test_view.dev  test_view.dev      [.] cram_to_bam
     5.93%           266  test_view.dev  libc.so.6          [.] __memmove_avx512_unaligned_erms
     3.13%           129  test_view.dev  libc.so.6          [.] __strncpy_evex
     2.50%            92  test_view.dev  [unknown]          [k] 0xffffffff82cb9b7e
     2.37%            97  test_view.dev  test_view.dev      [.] cram_get_seq
     2.25%            92  test_view.dev  test_view.dev      [.] sam_read1
     2.11%            92  test_view.dev  test_view.dev      [.] bam_cigar2rqlens
     1.65%            64  test_view.dev  test_view.dev      [.] cram_get_bam_seq

With this PR:

@ qpg-gpu-01[samtools.../htslib]; ./test/test_view.cram2bam -B -@32  ~/scratch/data/novaseq.10m.cram & perf record -t $! 2>/dev/null;perf report -n 2>/dev/null| head -20 | egrep -v '^#'
[1] 3634842
[1]+  Done                    ./test/test_view.cram2bam -B -@32 ~/scratch/data/novaseq.10m.cram
    15.65%           101  test_view.cram2  test_view.cram2bam  [.] cram_get_bam_seq
    11.31%            65  test_view.cram2  libc.so.6           [.] _int_free
     6.69%            38  test_view.cram2  [unknown]           [k] 0xffffffff829076fb
     5.17%            34  test_view.cram2  test_view.cram2bam  [.] cram_get_seq
     5.12%            29  test_view.cram2  [unknown]           [k] 0xffffffff8291f048
     4.73%            32  test_view.cram2  test_view.cram2bam  [.] sam_read1
     4.24%            24  test_view.cram2  [unknown]           [k] 0xffffffff828d0d41
     4.23%            24  test_view.cram2  [unknown]           [k] 0xffffffff8292a32a
     3.35%            19  test_view.cram2  [unknown]           [k] 0xffffffff82929c09

Wall clock benchmarks with 32 cores show a doubling in throughput (making it slightly faster than BAM, but that's due to the same main-thread excessive work problem there too):

$ for i in `seq 1 10`;do /usr/bin/time -f '%e %U %S'  ./test/test_view.dev 
-B -@32  ~/scratch/data/novaseq.10m.cram;done
1.14 10.34 0.76
1.06 9.13 0.69
1.11 9.42 0.72
1.12 10.20 0.62
1.00 8.72 0.57
1.23 8.85 0.77
1.46 9.13 0.79
1.17 8.67 0.63
1.30 8.73 0.69
0.99 9.64 0.54


$ for i in `seq 1 10`;do /usr/bin/time -f '%e %U %S'  ./test/test_view.cram2bam -B -@32  ~/scratch/data/novaseq.10m.cram;done
0.50 11.83 0.98
0.50 10.32 0.82
0.45 10.13 0.82
0.42 10.10 0.78
0.50 10.19 0.87
0.47 10.07 0.81
0.48 10.15 0.98
0.49 10.03 1.00
0.51 10.06 0.85
0.49 10.28 0.76

[daviesrob]

Looks like this does a good job of moving more work into the thread pool. Would it be worth doing the CG processing there too? I think it would just need bam_tag2cigar() to be made non-static so that it can be called from cram/cram_decode.c.

[jkbonfield]

Looks like this does a good job of moving more work into the thread pool. Would it be worth doing the CG processing there too? I think it would just need bam_tag2cigar() to be made non-static so that it can be called from cram/cram_decode.c.

Possibly, but I doubt this code is called very often. It's only an issue on ultra-long ONT reads, and maybe not even there it's not that often any more due to improved quality. If we moved it, we'd still need a flag to say whether or not it's been done already, so that's changed.

I think I'd need to find some test data to evaluate. (Most of mine got lost on lustre.) It may be an unnecessary optimisation though.

What's more important appears to be careful malloc handling. As we add more threads, we add more bam-list allocations. These appear to cause malloc heap thrashing. GNU malloc is particularly badly performing here. TCMalloc is double the speed at high thread count as it has minimal slow down with more threads. GNU malloc performs badly as it gives back memory and then has to reallocate it again. Defining MALLOC_MMAP_MAX_=0 MALLOC_TOP_PAD_=100000000 solves this, but in theory it's not meant to be giving the memory back anyway as this is what the bam-list bit was meant to solve. Maybe the memory thrashing is somewhere else. I need to diagnose that.

(Io_lib just used mallopt to work around gnu malloc performance, but it's not ideal in a library)

@ [samtools.../htslib]; for i in `seq 1 5`;do /usr/bin/time -f '%e %U %S'  ./test/test_view -B -@64
  ~/scratch/data/novaseq.10m.31.cram;done
3.27 28.82 6.30
3.17 28.65 6.35
3.30 28.32 6.44
3.32 28.42 6.67
3.27 29.00 6.53

@ [samtools.../htslib]; for i in `seq 1 5`;do MALLOC_MMAP_MAX_=0 MALLOC_TOP_PAD_=100000000 /usr/bin/time -f '%e %U %S'  ./test/test_view -B -@64  ~/scratch/data/novaseq.10m.31.cram;done
1.47 28.74 2.88
1.49 29.27 2.91
1.46 28.18 2.71
1.48 28.72 2.71
1.48 28.13 2.79

@ [samtools.../htslib]; for i in `seq 1 5`;do LD_PRELOAD=~/scratch/opt/lib/libtcmalloc.so /usr/bin/time -f '%e %U %S'  ./test/test_view -B -@64  ~/scratch/data/novaseq.10m.31.cram;done
1.66 35.33 3.25
1.61 34.49 3.24
1.61 34.72 3.07
1.64 34.43 3.29
1.65 34.65 3.35

(Note this machine has 16 cores each with 2 hyperthreads, so this is just a stress test of memory by over-committing on bam-list arrays.)

Edit: It's not the bam_list allocations, but the cram record one:

    if (!(s->crecs = hts_malloc_p(sizeof(*s->crecs), s->hdr->num_records)))     
        return -1;                                                              

This is probably another PR, but I'll investigate

[jkbonfield]

I have ideas for reducing malloc pressure, but it may come later depending on how long this sits in review (as it'll need to be based on top of these commits).

I've fixed the reallocs and removed needless code. I've rebased it too. I'll take a look at the bam_tag2cigar change as it may make the interface easier so perhaps it's worth it for that reason alone.

[jkbonfield]

I can squash when you wish, but figured it's helpful for review to not do that just now so you can see the new commits first. (Or squash it yourself if you wish.)

The bam_tag2cigar was quite simple to do.

[daviesrob]

Looks OK, apart from a couple of potential leaks on realloc().

[jkbonfield]

I went through all the reallocs in cram_encode.c and cram_decode.c and found a few more, so I fixed those too.
Now pushed up.

I also have another branch that creates a cram_list struct to go along the bam_list struct and reuses it between cram slices.
However it doesn't really seem to help much and we still need mallopt envs to coerce GNU malloc to not be so rubbish! (Or use tcmalloc). Frustrating. So for now I'm not adding that extra complexity.

I'm also now finding it hard to find any actual speed benefit on our modern machines (but it is faster on the old seq4 node). It's still provably faster by profiling the main thread, but something else is slowing it up now. Such are the vagaries of I/O and shared CPUs, but at least it can be faster when the conditions are favourable.