Feature request: bulk distributed access

[tomwhite]

Thank you for all your work on this library! We are using bgen-reader-py to implement a Dask-based BGEN reader in sgkit-dev/sgkit-bgen#1. There is already some support for Dask, but (as far as I can tell) the delayed object that is created is for a single row (variant), rather than a chunk of rows. (See https://github.com/limix/bgen-reader-py/blob/master/bgen_reader/_genotype.py#L46.) This is a good choice for random access, but may be prohibitive for bulk access.

Would it be possible to support bulk access with an API that returns the genotype probability array as a Dask array?

[CarlKCarlK]

Tom, Danilo (the lead on this project) and I have a beta for a 2nd API that supports bulk reads. It is as fast as we could make it (reading about 4 million distributions per second). It is based on NumPy rather than Dask, but I assume you could convert any data to Dask as needed. A beta quickstart document is here: https://fastlmm.github.io/PySnpTools/branch/bgen2/index.html. You can install the beta for Linux 3.7 from here (3.8 is where’d you expect. Windows and MacOS available, too.) pip install https://github.com/fastlmm/PySnpTools/releases/download/v4.0.4/bgen_reader-4.0.4-cp37-cp37m-manylinux1_x86_64.whl Let us know if this seems useful for you. * Carl https://github.com/carlkcarlk From: Tom White <notifications@github.com> Sent: Thursday, July 23, 2020 7:01 AM To: limix/bgen-reader-py <bgen-reader-py@noreply.github.com> Cc: Subscribed <subscribed@noreply.github.com> Subject: [limix/bgen-reader-py] Feature request: bulk distributed access (#30) Thank you for all your work on this library! We are using bgen-reader-py to implement a Dask-based BGEN reader in sgkit-dev/sgkit-bgen#1<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fpystatgen%2Fsgkit-bgen%2Fpull%2F1&data=02%7C01%7C%7C6632ef428cb84168572c08d82f10d6db%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637311096648889003&sdata=4NOVsKG19gUsgPJFUvAucIb6rKBxA7EY1xvF%2Fhdzogk%3D&reserved=0>. There is already some support for Dask, but (as far as I can tell) the delayed object that is created is for a single row (variant), rather than a chunk of rows. (See https://github.com/limix/bgen-reader-py/blob/master/bgen_reader/_genotype.py#L46<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Flimix%2Fbgen-reader-py%2Fblob%2Fmaster%2Fbgen_reader%2F_genotype.py%23L46&data=02%7C01%7C%7C6632ef428cb84168572c08d82f10d6db%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637311096648898996&sdata=j7wcgygIfgRrUUSdUzPoK5pZbdsINUjiHrlCTVzCEkw%3D&reserved=0>.) This is a good choice for random access, but may be prohibitive for bulk access. Would it be possible to support bulk access with an API that returns the genotype probability array as a Dask array? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Flimix%2Fbgen-reader-py%2Fissues%2F30&data=02%7C01%7C%7C6632ef428cb84168572c08d82f10d6db%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637311096648898996&sdata=zTNVSdCWbwfmRzy5AgyWeYEYl%2BI2UCJqGd5rPEhqdxw%3D&reserved=0>, or unsubscribe<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABR65P5NIAVH5B7ISVNJAOTR5A7B7ANCNFSM4PFYCRGQ&data=02%7C01%7C%7C6632ef428cb84168572c08d82f10d6db%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637311096648908993&sdata=GeYp1%2FTheXaZNvd2siDnt7JVvHSWF%2B46bx2i2%2BZ5210%3D&reserved=0>.

[tomwhite]

@CarlKCarlK that's great to know! That looks very useful. I will dig in and try it out soon.

[horta]

That should be done by end of this month (August).

[horta]

@tomwhite I'm trying to figure out a nice way to provide bulk access.

As you probably know, the probabilities matrix of genotype 0 maybe have different number of columns of genotype 1 as the ploidy is allowed to vary. Not only that but a genotype is (in bgen terms) defined in terms of probabilities as well as being phased, ploidy, missingness.

In [6]: genotype[0].compute()
Out[6]:
{'probs': array([[1., 0., 0.],
        [1., 0., 0.],
        [1., 0., 0.],
        ...,
        [1., 0., 0.],
        [1., 0., 0.],
        [1., 0., 0.]]),
 'phased': False,
 'ploidy': array([2, 2, 2, ..., 2, 2, 2]),
 'missing': array([False, False, False, ..., False, False, False])}

In [7]: genotype[1].compute()
Out[7]:
{'probs': array([[nan, nan, nan],
        [nan, nan, nan],
        [nan, nan, nan],
        ...,
        [ 1.,  0.,  0.],
        [ 1.,  0.,  0.],
        [nan, nan, nan]]),
 'phased': False,
 'ploidy': array([2, 2, 2, ..., 2, 2, 2]),
 'missing': array([ True,  True,  True, ..., False, False,  True])}

Without trying to be fancy, a quick way to implement bulk access is by indexing genotype by its "bulk index" instead of its genotype index. Something like that:

In [6]: bulk[0].compute().genotype[0]
Out[6]:
{'probs': array([[1., 0., 0.],
        [1., 0., 0.],
        [1., 0., 0.],
        ...,
        [1., 0., 0.],
        [1., 0., 0.],
        [1., 0., 0.]]),
 'phased': False,
 'ploidy': array([2, 2, 2, ..., 2, 2, 2]),
 'missing': array([False, False, False, ..., False, False, False])}

In [7]: bulk[0].compute().genotype[1]
Out[7]:
{'probs': array([[nan, nan, nan],
        [nan, nan, nan],
        [nan, nan, nan],
        ...,
        [ 1.,  0.,  0.],
        [ 1.,  0.,  0.],
        [nan, nan, nan]]),
 'phased': False,
 'ploidy': array([2, 2, 2, ..., 2, 2, 2]),
 'missing': array([ True,  True,  True, ..., False, False,  True])}

Does it appeal you?

[horta]

@CarlKCarlK pointed out that he already have bulk access implemented (in the numpy-like interface, accessed via open_bgen). I'm also interested in having that in the dask interface: via read_bgen

[hammer]

Just want to ensure @eric-czech is in this conversation too as he and Carl have been in discussion separately.

[eric-czech]

Hey @CarlKCarlK, I tested this on a single UKB file and had some results to share.

I ran a script for pybgen, bgen_reader, and the new beta to compare times / resource usage on a single XY PAR file scan. The code and the times I recorded are at https://github.com/eric-czech/sgkit-dev/tree/master/io/bgen/ukb-test.

The file is only 4.6G and here are the times I got on a single core:

• PyBGEN: ~28 minutes (the .bgi file already exists)
• bgen_reader: I stopped it after an hour but it was on track to take ~4 hours
• bgen_reader2: ~34 minutes not including the metadata construction time

I did a little more profiling on the new beta in bgen_reader2_prof.ipynb, though this only loads 10k variants of 45,906 to speed things up. I found that in reading these 10k variants after deleting the metafile, I see resource utilization like this:

The metafile apparently takes about 2m 40s to create, which would be approximately 8% of the overall time. It seems to take a somewhat marginal amount of time compared to the scan, and the cpu/memory usage profiles are otherwise very stable.

Is this about what you would have expected? Our use case would be to traverse the file once and convert it to another format -- is there any relatively simple way to speed up a scan like that?

A few smaller questions:

• Have you ever used tqdm? Assuming the custom progress prints are coming from python and not c++, that could make things easier. I noticed a small bug with each print where they seem to be off on the end interval, i.e.

metadata -- time=0:02:41.61, step 3: part 45,000 of 45,906
Found 10000 variants
reading -- time=0:00:34.72, part 990 of 1,000 # --> Never makes it to 1,000 
reading -- time=0:00:37.07, part 990 of 1,000

• I was doing my scan by grabbing row slices in batches of 1000 variants but it looked like the code was actually pulling them in smaller batches of 10. Each "part N" message was incrementing by 10 from 0 up to 1000. Does it work that way regardless of the file size?

[hammer]

If the progress prints are coming from C++, https://github.com/p-ranav/indicators is an option.

[CarlKCarlK]

Eric, Thanks! At the bottom of this node is the inner loop of the Python Bgen2 reader. In summary, it reads one variant at a time with these steps: * Open that variant * If the buffer can be reused (usually the case), do nothing, otherwise allocate a new buffer (for ‘complex’ files) * Read the probability data * Copy results into output array * (If samples are also being selected, do that, too) * Close the variant. The Open, Read, and Close methods are @Danilo Horta<mailto:danilo.horta@pm.me>’s C/C++ code. I haven’t profiled a run (we should), but I assume most of the time is being spent in lib.bgen_genotype_read. * Carl p.s. I fixed the logging off-by-one errors and added a TODO to look at tqdm<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Ftqdm%2Ftqdm&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965685853&sdata=2Ck9mdIa37hAV6Z4ComRBKkqFiDfeGda805JpliK1ak%3D&reserved=0>.. From https://github.com/limix/bgen-reader-py/blob/bgen2memmap/bgen_reader/_bgen2.py with _log_in_place("reading", self._verbose) as updater: for out_index, vaddr0 in enumerate(vaddr): if out_index % vaddr_per_second == 0: updater("part {0:,} of {1:,}".format(out_index, len(vaddr))) genotype = lib.bgen_file_open_genotype(self._bgen._bgen_file, vaddr0) if return_probabilities: if ( prob_buffer is None or ncombinations[out_index] != prob_buffer.shape[-1] ): prob_buffer = np.full( # LATER could offer an option to memmap this buffer (self.nsamples, ncombinations[out_index]), np.nan, order="C", dtype="float64", ) lib.bgen_genotype_read( genotype, ffi.cast("double *", prob_buffer.ctypes.data) ) val[:, out_index, : ncombinations[out_index]] = ( prob_buffer if (samples_index is self._sample_range) else prob_buffer[samples_index, :] ) if return_missings: […] if return_ploidies: […] lib.bgen_genotype_close(genotype) From: Jeff Hammerbacher <notifications@github.com> Sent: Monday, August 10, 2020 12:39 PM To: limix/bgen-reader-py <bgen-reader-py@noreply.github.com> Cc: Carl Kadie <carlk@msn.com>; Mention <mention@noreply.github.com> Subject: Re: [limix/bgen-reader-py] Feature request: bulk distributed access (#30) If the progress prints are coming from C++, https://github.com/p-ranav/indicators<https://eur06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fp-ranav%2Findicators&data=02%7C01%7C%7C468ab92806744667e3bd08d83d64fbea%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326851211525643&sdata=Fthje7kaMsdQmxJJPDsdpAHW0EeAooKyc5rp1xhhM60%3D&reserved=0> is an option. — You are receiving this because you were mentioned. Reply to this email directly From: Eric Czech <notifications@github.com> Sent: Monday, August 10, 2020 12:08 PM To: limix/bgen-reader-py <bgen-reader-py@noreply.github.com> Cc: Carl Kadie <carlk@msn.com>; Mention <mention@noreply.github.com> Subject: Re: [limix/bgen-reader-py] Feature request: bulk distributed access (#30) Hey @CarlKCarlK<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FCarlKCarlK&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965645877&sdata=sqBaai68dP2FhoJn%2BuedXbRih3MaXzElAo%2Fs35LpbR0%3D&reserved=0>, I tested this on a single UKB file and had some results to share. I ran a script for pybgen, bgen_reader, and the new beta to compare times / resource usage on a single XY PAR file scan. The code and the times I recorded are at https://github.com/eric-czech/sgkit-dev/tree/master/io/bgen/ukb-test<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Feric-czech%2Fsgkit-dev%2Ftree%2Fmaster%2Fio%2Fbgen%2Fukb-test&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965655878&sdata=KFJnv0cNdP3q9d2OreRZD29XLvyewpFz%2FOYWTxF%2BnOg%3D&reserved=0>. The file is only 4.6G and here are the times I got on a single core: * PyBGEN: ~28 minutes (the .bgi file already exists) * bgen_reader: I stopped it after an hour but it was on track to take ~4 hours * bgen_reader2: ~34 minutes not including the metadata construction time I did a little more profiling on the new beta in bgen_reader2_prof.ipynb<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fnbviewer.jupyter.org%2Fgithub%2Feric-czech%2Fsgkit-dev%2Fblob%2F34dba75a7db899b046d779e7373c4eb2ebd258ca%2Fio%2Fbgen%2Fukb-test%2Fbgen_reader2_prof.ipynb&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965665867&sdata=QgVmI7N7mrtJ1RvViuNMj4R9w5I7IXELbByS2ouLI3Y%3D&reserved=0>. I found that in reading 10k variants of 45,906 after deleting the metafile, I see resource utilization like this: [prof]<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Feric-czech%2Fsgkit-dev%2Fraw%2Fmaster%2Fio%2Fbgen%2Fukb-test%2Fbgen_reader2_prof.png&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965675861&sdata=lakJ68sPqEUEoOhScoGiMPQbRgnpb65yW6i0guS2XRc%3D&reserved=0> The metafile apparently takes about 2m 40s to create, which would be approximately 8% of the overall time. It seems to take a somewhat marginal amount of time compared to the scan, and the cpu/memory usage profiles are otherwise very stable. Is this about what you would have expected? Our use case would be to traverse the file once and convert it to another format -- is there any relatively simple way to speed up a scan like that? A few smaller questions: * Have you ever used tqdm<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Ftqdm%2Ftqdm&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965685853&sdata=2Ck9mdIa37hAV6Z4ComRBKkqFiDfeGda805JpliK1ak%3D&reserved=0>? Assuming the custom progress prints are coming from python and not c++, that could make things easier. I noticed a small bug with each print where they seem to be off on the end interval, i.e. metadata -- time=0:02:41.61, step 3: part 45,000 of 45,906 Found 10000 variants reading -- time=0:00:34.72, part 990 of 1,000 # --> Never makes it to 1,000 reading -- time=0:00:37.07, part 990 of 1,000 * I was doing my scan by grabbing row slices in batches of 1000 variants but it looked like the code was actually pulling them in smaller batches of 10. Each "part N" message was incrementing by 10 from 0 up to 1000. Does it work that way regardless of the file size? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Flimix%2Fbgen-reader-py%2Fissues%2F30%23issuecomment-671534711&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965695848&sdata=SsSP0Qz5tGe8ZoR3Yiar5D3a5YXAlJoES7t4SFZCMp4%3D&reserved=0>, or unsubscribe<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABR65P2MKYTB4YL64D2PG2TSABAR5ANCNFSM4PFYCRGQ&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965705845&sdata=dUOMyAnSB%2B8snkAI%2BpS87REdw%2FYq6ouDs5%2F3MN1%2F1zs%3D&reserved=0>.

[CarlKCarlK]

Update: Danilo has ideas on making it faster with C/C++ batch readers. I did the profiling: 84% of time was in the current read, 12% in the open. Everything else was a small %. * Carl p.s. Notes to myself: * python -m cProfile -o test.profile test_bgen2.py * snakeviz test.profile From: Carl KADIE Sent: Monday, August 10, 2020 1:31 PM To: limix/bgen-reader-py <reply@reply.github.com>; limix/bgen-reader-py <bgen-reader-py@noreply.github.com>; Danilo Horta <danilo.horta@pm.me> Cc: Mention <mention@noreply.github.com> Subject: RE: [limix/bgen-reader-py] Feature request: bulk distributed access (#30) Eric, Thanks! At the bottom of this node is the inner loop of the Python Bgen2 reader. In summary, it reads one variant at a time with these steps: * Open that variant * If the buffer can be reused (usually the case), do nothing, otherwise allocate a new buffer (for ‘complex’ files) * Read the probability data * Copy results into output array * (If samples are also being selected, do that, too) * Close the variant. The Open, Read, and Close methods are @Danilo Horta<mailto:danilo.horta@pm.me>’s C/C++ code. I haven’t profiled a run (we should), but I assume most of the time is being spent in lib.bgen_genotype_read. * Carl p.s. I fixed the logging off-by-one errors and added a TODO to look at tqdm<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Ftqdm%2Ftqdm&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965685853&sdata=2Ck9mdIa37hAV6Z4ComRBKkqFiDfeGda805JpliK1ak%3D&reserved=0>.. From https://github.com/limix/bgen-reader-py/blob/bgen2memmap/bgen_reader/_bgen2.py with _log_in_place("reading", self._verbose) as updater: for out_index, vaddr0 in enumerate(vaddr): if out_index % vaddr_per_second == 0: updater("part {0:,} of {1:,}".format(out_index, len(vaddr))) genotype = lib.bgen_file_open_genotype(self._bgen._bgen_file, vaddr0) if return_probabilities: if ( prob_buffer is None or ncombinations[out_index] != prob_buffer.shape[-1] ): prob_buffer = np.full( # LATER could offer an option to memmap this buffer (self.nsamples, ncombinations[out_index]), np.nan, order="C", dtype="float64", ) lib.bgen_genotype_read( genotype, ffi.cast("double *", prob_buffer.ctypes.data) ) val[:, out_index, : ncombinations[out_index]] = ( prob_buffer if (samples_index is self._sample_range) else prob_buffer[samples_index, :] ) if return_missings: […] if return_ploidies: […] lib.bgen_genotype_close(genotype) From: Jeff Hammerbacher <notifications@github.com<mailto:notifications@github.com>> Sent: Monday, August 10, 2020 12:39 PM To: limix/bgen-reader-py <bgen-reader-py@noreply.github.com<mailto:bgen-reader-py@noreply.github.com>> Cc: Carl Kadie <carlk@msn.com<mailto:carlk@msn.com>>; Mention <mention@noreply.github.com<mailto:mention@noreply.github.com>> Subject: Re: [limix/bgen-reader-py] Feature request: bulk distributed access (#30) If the progress prints are coming from C++, https://github.com/p-ranav/indicators<https://eur06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fp-ranav%2Findicators&data=02%7C01%7C%7C468ab92806744667e3bd08d83d64fbea%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326851211525643&sdata=Fthje7kaMsdQmxJJPDsdpAHW0EeAooKyc5rp1xhhM60%3D&reserved=0> is an option. — You are receiving this because you were mentioned. Reply to this email directly From: Eric Czech <notifications@github.com<mailto:notifications@github.com>> Sent: Monday, August 10, 2020 12:08 PM To: limix/bgen-reader-py <bgen-reader-py@noreply.github.com<mailto:bgen-reader-py@noreply.github.com>> Cc: Carl Kadie <carlk@msn.com<mailto:carlk@msn.com>>; Mention <mention@noreply.github.com<mailto:mention@noreply.github.com>> Subject: Re: [limix/bgen-reader-py] Feature request: bulk distributed access (#30) Hey @CarlKCarlK<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FCarlKCarlK&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965645877&sdata=sqBaai68dP2FhoJn%2BuedXbRih3MaXzElAo%2Fs35LpbR0%3D&reserved=0>, I tested this on a single UKB file and had some results to share. I ran a script for pybgen, bgen_reader, and the new beta to compare times / resource usage on a single XY PAR file scan. The code and the times I recorded are at https://github.com/eric-czech/sgkit-dev/tree/master/io/bgen/ukb-test<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Feric-czech%2Fsgkit-dev%2Ftree%2Fmaster%2Fio%2Fbgen%2Fukb-test&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965655878&sdata=KFJnv0cNdP3q9d2OreRZD29XLvyewpFz%2FOYWTxF%2BnOg%3D&reserved=0>. The file is only 4.6G and here are the times I got on a single core: * PyBGEN: ~28 minutes (the .bgi file already exists) * bgen_reader: I stopped it after an hour but it was on track to take ~4 hours * bgen_reader2: ~34 minutes not including the metadata construction time I did a little more profiling on the new beta in bgen_reader2_prof.ipynb<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fnbviewer.jupyter.org%2Fgithub%2Feric-czech%2Fsgkit-dev%2Fblob%2F34dba75a7db899b046d779e7373c4eb2ebd258ca%2Fio%2Fbgen%2Fukb-test%2Fbgen_reader2_prof.ipynb&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965665867&sdata=QgVmI7N7mrtJ1RvViuNMj4R9w5I7IXELbByS2ouLI3Y%3D&reserved=0>. I found that in reading 10k variants of 45,906 after deleting the metafile, I see resource utilization like this: [prof]<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Feric-czech%2Fsgkit-dev%2Fraw%2Fmaster%2Fio%2Fbgen%2Fukb-test%2Fbgen_reader2_prof.png&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965675861&sdata=lakJ68sPqEUEoOhScoGiMPQbRgnpb65yW6i0guS2XRc%3D&reserved=0> The metafile apparently takes about 2m 40s to create, which would be approximately 8% of the overall time. It seems to take a somewhat marginal amount of time compared to the scan, and the cpu/memory usage profiles are otherwise very stable. Is this about what you would have expected? Our use case would be to traverse the file once and convert it to another format -- is there any relatively simple way to speed up a scan like that? A few smaller questions: * Have you ever used tqdm<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Ftqdm%2Ftqdm&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965685853&sdata=2Ck9mdIa37hAV6Z4ComRBKkqFiDfeGda805JpliK1ak%3D&reserved=0>? Assuming the custom progress prints are coming from python and not c++, that could make things easier. I noticed a small bug with each print where they seem to be off on the end interval, i.e. metadata -- time=0:02:41.61, step 3: part 45,000 of 45,906 Found 10000 variants reading -- time=0:00:34.72, part 990 of 1,000 # --> Never makes it to 1,000 reading -- time=0:00:37.07, part 990 of 1,000 * I was doing my scan by grabbing row slices in batches of 1000 variants but it looked like the code was actually pulling them in smaller batches of 10. Each "part N" message was incrementing by 10 from 0 up to 1000. Does it work that way regardless of the file size? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Flimix%2Fbgen-reader-py%2Fissues%2F30%23issuecomment-671534711&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965695848&sdata=SsSP0Qz5tGe8ZoR3Yiar5D3a5YXAlJoES7t4SFZCMp4%3D&reserved=0>, or unsubscribe<https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABR65P2MKYTB4YL64D2PG2TSABAR5ANCNFSM4PFYCRGQ&data=02%7C01%7C%7C1aa27b04dcc5462cbb6508d83d60bc54%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326832965705845&sdata=dUOMyAnSB%2B8snkAI%2BpS87REdw%2FYq6ouDs5%2F3MN1%2F1zs%3D&reserved=0>.

[horta]

Thanks for the benchmark, @eric-czech !

I have found some bootlenecks on the bgen-reader. Here is a small comparison now:

% python pybgen_exec.py run --path=$HOME/carl/merged_487400x220000.bgen
  0%|          | 778/220000 [00:18<1:28:18, 41.37it/s]

% python bgen_reader_exec.py run --path=$HOME/carl/merged_487400x220000.bgen
  0%|          | 483/220000 [00:10<1:22:19, 44.44it/s]

I have stopped more or less after 10 seconds just to have an idea.

I modified your pybgen_exec.py only to use tqdm:

import fire
import tqdm
from pybgen import PyBGEN


def run(path):
    ct = 0
    with PyBGEN(path) as bgen:
        n = bgen.nb_variants
        for info, dosage in tqdm.tqdm(bgen, total=n):
            assert dosage.ndim == 1
            ct += dosage.size
    print(f"Number of entries read: {ct}")


fire.Fire()

I will try to release tonight a version 4.0.5 for bgen-reader with that speed-up. It would be great if you could find time to run the benchmarks again on that new version (I come back here as soon as I release it).

It will come with the Carl API as well, which might be in the state of experimental but I think that is fine for now.

Notice that I've not used bulk reading (nor have implemented it yet), although reading the variants in order using bgen_reader.read_bgen is faster than in random order.

[horta]

@hammer Thanks for the pointer. I didnt know about that library.

I've coded the low-level of bgen in C and at that time I did not find a nice progress bar library for C only. So I coded mine: https://github.com/horta/almosthere
It is very humble compared to the one you showed: ASCII based only.

[CarlKCarlK]

Aside: I have my own, too 😊 A single 45-lines Python function called “_log_in_place". Here is a usage example. Notice the “i+1”. That fixes the off-by-one error Eric found. with _log_in_place("metadata", self._verbose) as updater: for i in range(self.nsamples): if i % 1000 == 0: updater("'sample range': part {0:,} of {1:,}".format(i+1, self.nsamples)) sample_range_memmap[i] = i * Carl From: Danilo Horta <notifications@github.com> Sent: Monday, August 10, 2020 3:02 PM To: limix/bgen-reader-py <bgen-reader-py@noreply.github.com> Cc: Carl Kadie <carlk@msn.com>; Mention <mention@noreply.github.com> Subject: Re: [limix/bgen-reader-py] Feature request: bulk distributed access (#30) @hammer<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fhammer&data=02%7C01%7C%7C74de41b3873447dfb6f608d83d79075c%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326937302180578&sdata=WZVbsv7QfG6QSjZFnVwEgJMoXwGVOEeRYUZP4V2oS1I%3D&reserved=0> Thanks for the pointer. I didnt know about that library. I've coded the low-level of bgen in C and at that time I did not find a nice progress bar library for C only. So I coded mine: https://github.com/horta/almosthere<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fhorta%2Falmosthere&data=02%7C01%7C%7C74de41b3873447dfb6f608d83d79075c%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326937302180578&sdata=J9YQZijOG%2F%2FXAAtc1f9TEMZatyYWBvrUDloorhioLBo%3D&reserved=0> It is very humble compared to the one you showed: ASCII based only. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Flimix%2Fbgen-reader-py%2Fissues%2F30%23issuecomment-671611940&data=02%7C01%7C%7C74de41b3873447dfb6f608d83d79075c%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326937302190564&sdata=G0FrivExKUsOtQ1DRMMmVfaVAkJYA9xucTSHF1kTTdc%3D&reserved=0>, or unsubscribe<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABR65PYV5BIHIW4FEOFXQSTSABU6BANCNFSM4PFYCRGQ&data=02%7C01%7C%7C74de41b3873447dfb6f608d83d79075c%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637326937302190564&sdata=qeFoVB0hyQg1kl0TQRHnHewAHfIZNSPAwbHn7ICpVlk%3D&reserved=0>.