My tsinfer inference is taking too long when matching ancestors

[stsmall]

Moving discussion from slack channel to here.
Issue: tsinfer is taking a long-time to run the ma step using data with many sites (5e6) and individuals (~1500).
some suggestions:
AncestorData.truncate_ancestors()
time-binning (#284) using https://github.com/awohns/unified_genealogy_paper/blob/master/all-data/bin_missing.py

[stsmall]

Following suggestions of time binning and truncating ancestors (un-dated tree, low=0, high=0.5) ... is faster from 56s/it down to 26s/it. Still not making much progress though.
ma-match (1/1) 0%|▎ | 11.8k/2.54M [11:27:43, 26.7s/it, epoch=0.9989775051124744, nanc=4264]

[stsmall]

I discovered that the bin_missing.py does not do any binning when there is no missing data. I altered the script and reduced the times.

bin_missing.py ::  Number of samples: 2934 . Number of discrete times: 2933
bin_freq.py (rewrite) :: Number of samples: 2934 . Number of discrete times: 1001

[jeromekelleher]

Would you mind pasting in more details here please @stsmall? The output of python -m tsinfer ls <file> for your sample data and ancestor data files would be a good start. Also if you could show the output of your ma process when you have verbosity turned up a bit would help too.

[jeromekelleher]

Ah, I see. You can get the logging to emit using something like this in your script:

import daiquiri

daiquiri.setup(level="INFO")

you can set level to DEBUG also, but it's a lot of info and probably not helpful.

[stsmall]

I was able to run the CLI to completion (see below).
However, we were hoping to use a recombination rate or map rate in the analyses. I am currently running that analyses, ma, with the logging turned on. It has made about 1% progress in the last 20 hours.

ma-match (1/1)100%|| 2.54M/2.54M [26:35:23, 26.5it/s, epoch=0.0006816632583503749, nanc=633422]
ms-match (1/3)100%|| 2.93k/2.93k [7:55:32, 9.72s/it]
ms-paths (2/3)100%|| 2.93k/2.93k [02:39, 18.4it/s]
ms-muts  (3/3)100%|| 2.76M/2.76M [00:51, 53.6kit/s]
ms-xsites (4/3)100%|| 2.52M/2.52M [3:28:23, 202it/s]

path      = gamb.X-bin.ma.tree
size      = 1.8 GiB
edges     = 45350685
trees     = 2246244
sites     = 2757046
mutations = 2757046

path      = gamb.X-bin.ms.tree
size      = 1.6 GiB
edges     = 34220084
trees     = 2441577
sites     = 5278617
mutations = 5278617

[jeromekelleher]

It's debateable whether there's any real benefit to running match-ancestors and non-zero mutation rates. It's not obvious what the interpretation is, given that we're running with synthetically generated ancestors.

I think a good approach here would be to subset your data down to say, 5% of a "representative" chunk of the chromosome, and then do some exploratory work on the mismatch parameters. Then once you've figured out what combination of parameters seems best, try to scale that up to the whole chr.

My prior would be that having a mismatch rate of 0 for match-ancestors is a perfectly reasonably prior, though.

[stsmall]

Thanks for the advice Jerome!
Just so I understand.
you are suggesting to run match-ancestors with no recombination and no mismatch set and then add the recombination rate and mismatch for the match-samples.

Are you then suggesting I do this for the representative chunks (same workflow) to find the appropriate mismatch for the match-samples or to add into the match-ancestors?

Also you mention non-zero mutation rate ... do you mean recombination rate or am I missing something?

[jeromekelleher]

I'd start with a smaller dataset in order to figure out what the appropriate parameters are. So, say 10% of the variants on one chromosome.

I don't think specifying a varying recombination rate along the genome is going to have much effect on match-ancestors, as the breakpoints are usually driven by the end-points of haplotypes (I could well be wrong about this, though). Similarly, I would think that allowing mismatch during match-ancestors is of limited benefit, so setting mismatch=0 for that phase should be fine.

Adding nonzero mismatch and varying recombination rates for match-samples is much more useful, I think, and much faster because there's relatively many fewer samples than ancestors.

Does that clear things up a bit?

[stsmall]

yes, np. Original samples file.

path                  = gamb.X.samples
file_size             = 178.2 MiB
format_name           = tsinfer-sample-data
format_version        = (5, 0)
finalised             = True
uuid                  = 9f7d4c4d-5fea-4f9b-a892-e07ca1610d17
num_provenances       = 0
provenances/timestamp = shape=(0,); dtype=object;
provenances/record    = shape=(0,); dtype=object;
sequence_length       = 24393108.0
metadata_schema       = {'additionalProperties': True, 'codec': 'json', 'properties': {}, 'required': [], 'type': 'object'}
metadata              = {}
num_populations       = 100
num_individuals       = 1467
num_samples           = 2934
num_sites             = 5278617
populations/metadata_schema = None
populations/metadata  = shape=(100,); dtype=object;
individuals/metadata_schema = None
individuals/metadata  = shape=(1467,); dtype=object;
individuals/location  = shape=(1467,); dtype=object;
individuals/time      = shape=(1467,); dtype=float64;uncompressed size=11.7 kB
individuals/population = shape=(1467,); dtype=int32;uncompressed size=5.9 kB
individuals/flags     = shape=(1467,); dtype=uint32;uncompressed size=5.9 kB
samples/individual    = shape=(2934,); dtype=int32;uncompressed size=11.7 kB
sites/position        = shape=(5278617,); dtype=float64;uncompressed size=42.2 MB
sites/time            = shape=(5278617,); dtype=float64;uncompressed size=42.2 MB
sites/alleles         = shape=(5278617,); dtype=object;
sites/genotypes       = shape=(5278617, 2934); dtype=int8;uncompressed size=15.5 GB
sites/metadata_schema = None
sites/metadata        = shape=(5278617,); dtype=object;

[stsmall]

after running bin_missing.py. Also the ouput message of bin_missing.py
Number of samples: 2934 . Number of discrete times: 2933

path                  = gamb.X-bin.samples
file_size             = 189.8 MiB
format_name           = tsinfer-sample-data
format_version        = (5, 0)
finalised             = True
uuid                  = d2b83bc6-1d7f-4259-9188-ace0ba8c2280
num_provenances       = 0
provenances/timestamp = shape=(0,); dtype=object;
provenances/record    = shape=(0,); dtype=object;
sequence_length       = 24393108.0
metadata_schema       = {'additionalProperties': True, 'codec': 'json', 'properties': {}, 'required': [], 'type': 'object'}
metadata              = {}
num_populations       = 100
num_individuals       = 1467
num_samples           = 2934
num_sites             = 5278617
populations/metadata_schema = None
populations/metadata  = shape=(100,); dtype=object;
individuals/metadata_schema = None
individuals/metadata  = shape=(1467,); dtype=object;
individuals/location  = shape=(1467,); dtype=object;
individuals/time      = shape=(1467,); dtype=float64;uncompressed size=11.7 kB
individuals/population = shape=(1467,); dtype=int32;uncompressed size=5.9 kB
individuals/flags     = shape=(1467,); dtype=uint32;uncompressed size=5.9 kB
samples/individual    = shape=(2934,); dtype=int32;uncompressed size=11.7 kB
sites/position        = shape=(5278617,); dtype=float64;uncompressed size=42.2 MB
sites/time            = shape=(5278617,); dtype=float64;uncompressed size=42.2 MB
sites/alleles         = shape=(5278617,); dtype=object;
sites/genotypes       = shape=(5278617, 2934); dtype=int8;uncompressed size=15.5 GB
sites/metadata_schema = None
sites/metadata        = shape=(5278617,); dtype=object;

[stsmall]

after running ga, then truncate_ancestors(low=0.0, high=0.5). These ranges were arbitrary as I was just testing runtimes.

path                  = gamb.X-bin.trees.ancestors
file_size             = 2.7 GiB
format_name           = tsinfer-ancestor-data
format_version        = (3, 0)
finalised             = True
uuid                  = 8bc7934d-44b0-4882-9f48-5f123e0118cc
num_provenances       = 1
provenances/timestamp = shape=(1,); dtype=object;
provenances/record    = shape=(1,); dtype=object;
sequence_length       = 24393108.0
sample_data_uuid      = d2b83bc6-1d7f-4259-9188-ace0ba8c2280
num_ancestors         = 2535144
num_sites             = 2757046
sites/position        = shape=(2757046,); dtype=float64;uncompressed size=22.1 MB
ancestors/start       = shape=(2535144,); dtype=int32;uncompressed size=10.1 MB
ancestors/end         = shape=(2535144,); dtype=int32;uncompressed size=10.1 MB
ancestors/time        = shape=(2535144,); dtype=float64;uncompressed size=20.3 MB
ancestors/focal_sites = shape=(2535144,); dtype=object;
ancestors/haplotype   = shape=(2535144,); dtype=object;