Trying to speed up STIR acquisition data algebra by avoiding unnecessary 0-fill of newly created data.

[evgueni-ovtchinnikov]

Changes in this pull request

Unnecessary 0-fill of newly created acquisition data reduced.

Testing performed

The following SIRF Python script (employing STIR via SIRF Python interface) was used for comparing STIR's ProjDataInMemory algebra performance with numpy:

import sirf.STIR as pet
pet.AcquisitionData.set_storage_scheme('memory')

acq = pet.AcquisitionData('prompts.hs')
x = acq.clone()
y = x*2

start = timeit.default_timer()
z = x + y
elapsed = timeit.default_timer() - start
print(f'z = x + y time {elapsed}')

start = timeit.default_timer()
x_arr = x.as_array()
y_arr = y.as_array()
elapsed = timeit.default_timer() - start
print(f'as_array time {elapsed}')
numpy_time = elapsed

start = timeit.default_timer()
z_arr = x_arr + y_arr
elapsed = timeit.default_timer() - start
print(f'numpy z = x + y time {elapsed}')
numpy_time += elapsed

start = timeit.default_timer()
z.fill(z_arr)
elapsed = timeit.default_timer() - start
print(f'fill time {elapsed}')
numpy_time += elapsed
print(f'total numpy z = x + y time {numpy_time}')

With a recent STIR, the output looks like

z = x + y time 0.5256878180002786
as_array time 0.23524911599997722
numpy z = x + y time 0.1330225659999087
fill time 0.07107986500022889
total numpy z = x + y time 0.4393515470001148

demonstrating that STIR's ProjDataInMemory algebra is much less efficient than numpy algebra applied to ProjDataInMemory data copied into numpy array. The investigation of this issue revealed that a considerable chunk of the computation time was wasted on zero-filling the newly created ProjDataInMemory object where x + y data is to end up.

This PR succeeded in reducing such waste, so that the timings have become

z = x + y time 0.24671401699015405
as_array time 0.23049703100696206
numpy z = x + y time 0.13044680299935862
fill time 0.0659959859913215
total numpy z = x + y time 0.4269398199976422

Related issues

Checklist before requesting a review

• [] I have performed a self-review of my code
• [] I have added docstrings/doxygen in line with the guidance in the developer guide
• [] I have implemented unit tests that cover any new or modified functionality (if applicable)
• The code builds and runs on my machine
• [] documentation/release_XXX.md has been updated with any functionality change (if applicable)

[KrisThielemans]

Ubuntu jobs fail, see #1550

MacOS is unaffected by that, but has

ERROR: NumericVectorWithOffset::xapyb: index ranges don't match

(note that this will be seen in Debug mode only. Otherwise it will likely lead to a segfault)

[KrisThielemans]

While this would work, it could be a bit surprising to the user, as the effect of set_initialise_with_zeros() would be permanent. In the use-case that you have in mind, that'd be fine, but it'd have to come with serious warnings in the doxygen.

I'm rather in favour of having an explicit bool argument to grow and resize (defaulting to true). We could then have that for the constructor as well.

[KrisThielemans]

Please also use pre-commit, see https://github.com/UCL/STIR/blob/master/documentation/devel/git-hooks.md

[KrisThielemans]

No point in pushing until I've sorted out #1552, which you'll have to merge here.

[KrisThielemans]

#1552 is now merged, so merge master back on this branch.

[KrisThielemans]

I'm rather in favour of having an explicit bool argument to grow and resize (defaulting to true). We could then have that for the constructor as well.

I still think this would be a better option. @danieldeidda @markus-jehl what do you think?

[markus-jehl]

I tend to agree with @KrisThielemans, that explicitly stating it may be cleaner rather than relying on an internal setting on individual instances of the Array class. Just because I can imagine running into weird issues and struggling for a while before noticing what is going on 😅

[KrisThielemans]

I have some small comments, but realise now a conceptual difficulty: you had to add an arg to the virtual VectorWithOffset::grow. However, there it is unused (and would necessarily make sense, as "0" might not work for an arbitrary type). This is therefore still surprising to the user.

All this is more complicated due to the fact that adding default values to virtual functions is not possible and changing default values for virtual functions isn't recommended, see e.g. https://softwareengineering.stackexchange.com/questions/355074/is-overriding-a-pure-virtual-function-with-default-arguments-good-or-bad.

Could we do something like this

// make a "no_op" function object
virtual VectorWithOffset::grow(min_idx, max_idx, func& f); // apply `f` on new elements
virtual VectorWithOffset::grow(min_idx, max_idx); // no `f`, i.e. current version

virtual Array::grow(min_idx, max_idx); // assign with 0 (i.e. current version)
using virtual Array::grow(min_idx, max_idx, f); // use base-class version

if we don't want to initialise, then we could pass the no_op function.

Difficulties:

• it's a bit ugly to have to pass no_op (caused by the current Array::grow default being the wrong thing)
• to avoid overhead, we'd have to detect if f == no_op and avoids loops.

Not sure how to resolve this, unfortunately. Any ideas?

[KrisThielemans]

Suggested change

	inline virtual void grow(const int min_index, const int max_index, bool initialise_with_0 = true);
	inline void grow(const int min_index, const int max_index, bool initialise_with_0 = true) override;

[KrisThielemans]

Suggested change

	inline virtual void resize(const int min_index, const int max_index, bool initialise_with_0 = true);
	inline void resize(const int min_index, const int max_index, bool initialise_with_0 = true) override;

[KrisThielemans]

will need to be removed

[KrisThielemans]

will need to be removed