Computing Hessian products with image data objects and related stuff #1253

evgueni-ovtchinnikov · 2024-05-10T13:02:30Z

Changes in this pull request

Testing performed

Related issues

Fixes #1244, #1249

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
The code builds and runs on my machine
CHANGES.md has been updated with any functionality change

Contribution Notes

Please read and adhere to the contribution guidelines.

Please tick the following:

The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in SIRF (the Work) under the terms and conditions of the Apache-2.0 License.

…ted)

…sted)

KrisThielemans · 2024-05-10T13:38:50Z

src/xSTIR/cSTIR/cstir.cpp

+typedef xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF LMObjFun;
+
+extern "C"
+void* cSTIR_objFunListModeSetInterval(void* ptr_f, size_t ptr_data)


Not really for a question for this PR, but why size_t ptr_data and not void * ptr_data? On some systems, they are not the same size, and this could therefore create trouble.

I vaguely remember having some SWIG trouble with void* but no longer remember what it was. Perhaps with the latest SWIG void* is now ok.

BTW SWIG actually still does not like passing the numpy array data pointer as void* - tried with cSTIR_getImageData, got this:

Traceback (most recent call last): File "/home/sirfuser/devel/buildVM/sources/SIRF/examples/Python/PET/acquisition_data.py", line 172, in <module> main() File "/home/sirfuser/devel/buildVM/sources/SIRF/examples/Python/PET/acquisition_data.py", line 141, in main image_array = image.as_array() File "/home/sirfuser/devel/install/python/sirf/STIR.py", line 611, in as_array try_calling(pystir.cSTIR_getImageData(self.handle, array.ctypes.data)) File "/home/sirfuser/devel/install/python/sirf/pystir.py", line 306, in cSTIR_getImageData return _pystir.cSTIR_getImageData(ptr, ptr_data) TypeError: in method 'cSTIR_getImageData', argument 2 of type 'void *' sirfuser@vagrant:~/devel/buildVM/sources/SIRF/examples/Python/PET$

same story with float*:

Traceback (most recent call last): File "/home/sirfuser/devel/buildVM/sources/SIRF/examples/Python/PET/acquisition_data.py", line 172, in <module> main() File "/home/sirfuser/devel/buildVM/sources/SIRF/examples/Python/PET/acquisition_data.py", line 141, in main image_array = image.as_array() File "/home/sirfuser/devel/install/python/sirf/STIR.py", line 611, in as_array try_calling(pystir.cSTIR_getImageData(self.handle, array.ctypes.data)) File "/home/sirfuser/devel/install/python/sirf/pystir.py", line 306, in cSTIR_getImageData return _pystir.cSTIR_getImageData(ptr, ptr_data) TypeError: in method 'cSTIR_getImageData', argument 2 of type 'float *'

hmmm. weird. Here's an example that does that https://stackoverflow.com/a/37308401. Anyway, let's leave that for later!

src/xSTIR/cSTIR/cstir.cpp

KrisThielemans · 2024-05-10T13:44:06Z

src/xSTIR/cSTIR/cstir.cpp

+			fun.accumulate_sub_Hessian_times_input(output, curr_est, input, subset);
+		else {
+			int nsub = fun.get_num_subsets();
+			output.fill(0.0);


if we're accumulating in SIRF (which I don't think we should), then this fill is wrong. On the other hand, if we're not accumulating, then you'll need it also in the case subset>=0

moved the fill above if for now - will get rid of it when we switch to multiply_with_Hessian

KrisThielemans · 2024-05-10T13:46:17Z

Thanks, this already looks usable (with a caveat for the "accumulation"). Is this in a state that @gschramm can already try?

evgueni-ovtchinnikov · 2024-05-10T14:58:08Z

We still have #1251, which is due to the fact that PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin does not override empty ObjectiveFunction.get_subset_sensitivity(). I am afraid I cannot be of much help here.

… Python

gschramm · 2024-05-12T16:52:52Z

@KrisThielemans @evgueni-ovtchinnikov Just implemented a short test script for the PoissonLogL .accumulate_Hessian_times_input() here.

current_estimate = ref_recon.copy()
input_img = acq_data.create_uniform_image(value=1, xy=nxny)
np.random.seed(0)
input_img.fill(np.random.rand(*input_img.shape) * (obj_fun.get_subset_sensitivity(0).as_array() > 0) * current_estimate.max())

hess_out_img = acq_data.create_uniform_image(value=1.0, xy=nxny)

obj_fun.accumulate_Hessian_times_input(current_estimate, input_img, subset=0, out=hess_out_img)
hess_out_img2 = obj_fun.accumulate_Hessian_times_input(current_estimate, input_img, subset=0)

fig, ax = plt.subplots(1, 4, figsize=(16, 4), tight_layout=True)
ax[0].imshow(current_estimate.as_array()[71, :, :], cmap = 'Greys')
ax[1].imshow(input_img.as_array()[71, :, :], cmap = 'Greys')
ax[2].imshow(hess_out_img.as_array()[71, :, :], cmap = 'Greys')
ax[3].imshow(hess_out_img2.as_array()[71, :, :], cmap = 'Greys')
ax[0].set_title('current estimate', fontsize = 'medium')
ax[1].set_title('input', fontsize = 'medium')
ax[2].set_title('hess_out_img', fontsize = 'medium')
ax[3].set_title('hess_out_img2', fontsize = 'medium')
fig.show()

Unfortunately, the output of both calls (without and with using the out kwarg), don't seem to work.
hess_out_img is still an image full of ones and hess_out_img2 seems to be the same as input_img.

These results used:
STIR
commit feb6d85eadb392f5b8278d3b97ae2ee67ca439d9 (HEAD, origin/master, origin/HEAD, master)
Author: Kris Thielemans [email protected]
Date: Thu May 9 23:39:17 2024 +0100
SIRF
commit 66c35c7 (HEAD, origin/acc-hess)
Author: Evgueni Ovtchinnikov [email protected]
Date: Sat May 11 07:35:35 2024 +0000

KrisThielemans · 2024-05-12T19:32:01Z

This is a test with protection data, correct?

Also use auto more

gschramm · 2024-05-12T20:25:34Z

Indeed. The obj_fun is

obj_fun = sirf.STIR.make_Poisson_loglikelihood(acq_data)

KrisThielemans · 2024-05-12T20:32:16Z

@gschramm could you try again?

gschramm · 2024-05-12T20:53:33Z

Looks much better now. Tmr morning I will test whether we images are what we expect from the formulas.

KrisThielemans · 2024-05-12T21:09:26Z

Could work for listmode as well

gschramm · 2024-05-13T07:31:20Z

@KrisThielemans I just compared against manually computing the "Hessian multiply" calling the fwd / back projections in SIRF. The result is very close in the foreground (where the current estimate is >> 0), but not super close in the background (where the current estimate is close to 0). Any ideas why that is?

hess_out_img = obj_fun.accumulate_Hessian_times_input(current_estimate, input_img, subset=0)

# %%
# calcuate the the Hessian multiply "manually"

acq_model.set_up(acq_data, initial_image)
acq_model.num_subsets = num_subsets
acq_model.subset_num = 0

# get the linear (Ax) part of the Ax + b affine acq. model
lin_acq_model = acq_model.get_linear_acquisition_model()
lin_acq_model.num_subsets = num_subsets
lin_acq_model.subset_num = 0

# for the Hessian "multiply" we need the linear part of the acquisition model applied to the input image
input_img_fwd = lin_acq_model.forward(input_img)
current_estimate_fwd = acq_model.forward(current_estimate)
h = -acq_model.backward(acq_data*input_img_fwd / (current_estimate_fwd*current_estimate_fwd))

The Hessian outputs are shown using a "narrow" color window (top) and a wide color window (bottom) to show the differences in foreground and background.

KrisThielemans · 2024-05-13T08:38:13Z

hmmm. We try to cancel singularities in the division
https://github.com/UCL/STIR/blob/feb6d85eadb392f5b8278d3b97ae2ee67ca439d9/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx#L1128 with STIR's divide_and_truncate function. However, possibly that goes a bit wrong, as it was designed for (measured / estimated), while here it is used for (fwd * measured / estimated^2), so maybe thresholds are a bit off.

Of course, ideally there wouldn't be any singularities. Is your background strictly positive?

gschramm · 2024-05-13T09:52:14Z

The background term has indeed a few 0s (Siemens mMR data which I guess has virtual 0-sens. LORs).
In the meantime I also tested the listmode objective hessian. (i) currently, it returns the negative of what we want (output is positive, but should be negative). (ii) The "listmode hessian" is very close to my "manual sinogram" hessian (see below - I plotted the "minus the LM hessian output").

KrisThielemans · 2024-05-13T10:03:36Z

that's encouraging.

can you try your manual with some handling of 0/0 or so? A bit surprising it didn't generate NaNs then I guess. Easiest is to add an eps to the denominator.

I'll flip the sign of the LM Hessian. 😊

gschramm · 2024-05-13T10:47:21Z

adding an eps = 1e-8 in the denominator of the ratio to be back projected in the Hessian does not change much - still much closer to the LM Hessian applied.

KrisThielemans · 2024-05-13T11:47:08Z

ah well, this looks like a STIR issue then for the Hessian with projection data. Best to open an issue there.

KrisThielemans · 2024-05-13T11:48:09Z

@evgueni-ovtchinnikov can you then add multiply_with_Hessian? (Can just as well leave the accumulate version in I guess)

KrisThielemans · 2024-05-13T11:54:35Z

@gschramm what computation times are you getting?

evgueni-ovtchinnikov · 2024-05-13T11:59:19Z

multiply_with_Hessian: @KrisThielemans which version/branch of STIR should I use? I have rel_6.0.0 on my latest VMs.

gschramm · 2024-05-13T12:32:21Z

@gschramm what computation times are you getting?

Using 21 subsets, mMR LM file with 4,332,816 counts and acq_model = sirf.STIR.AcquisitionModelUsingRayTracingMatrix() and 1 tangential LOR:

In [5]: %timeit hess_out_img = obj_fun.accumulate_Hessian_times_input(current_estimate, input_img, subset=0)
2.7 s ± 45.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [6]: %timeit hess_out_img_lm = lm_obj_fun.accumulate_Hessian_times_input(current_estimate, input_img, subset=0)
9.43 s ± 249 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

KrisThielemans · 2024-05-13T12:34:48Z

multiply_with_Hessian: @KrisThielemans which version/branch of STIR should I use? I have rel_6.0.0 on my latest VMs.

you can use 6.0.0. Just add an extra member (that calls accumulate) and remove the fill in `accumulate.

KrisThielemans · 2024-05-13T12:39:48Z

Thanks @gschramm would you mind adding the timings in the STIR PR? Ideally also add gradient timings. I'd be interested as well in seeing what it does if num_subsets decreases (i.e. is the listmode calculations slowing us down, or the reading). Of course, you could also enable the lm_obj_fun caching, but this is going beyond what you care about. Thanks a lot for all your help.

gschramm · 2024-05-15T05:51:30Z

@KrisThielemans @evgueni-ovtchinnikov is accumulate_Hessian_times_input the final name? (I remember discussions on the naming) If possible, I'd like to finish my notebooks today.

KrisThielemans · 2024-05-15T06:47:15Z

@evgueni-ovtchinnikov is working on implementing #1244 (comment). Hopefully this will be done soon. @evgueni-ovtchinnikov if you have trouble, please ask for help.

KrisThielemans · 2024-05-15T06:49:41Z

@evgueni-ovtchinnikov here's the STIR test I talked about. https://github.com/UCL/STIR/blob/90cff236be35628f447ded4d98a046d5c4e3b316/src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h#L169
Best to implement the name first though.

see osem_reconstruction.py for a simple regression test for the Hessian multiplication

…(x)dx

KrisThielemans · 2024-05-15T16:59:36Z

@evgueni-ovtchinnikov you seem to be pretty close. Can you summarise current situation?

evgueni-ovtchinnikov · 2024-05-15T17:09:43Z

@KrisThielemans @gschramm I believe this PR can be merged - any objections?

evgueni-ovtchinnikov · 2024-05-15T17:12:16Z

builds failed at Coveralls step - any idea anyone why all of a sudden?

KrisThielemans

I think there are problems with the fill(0). It should NOT be done for the accumulate function, but has to be done for the multiply version.

Please batch commits of suggestions (I think via the Files tab), or do them manually of course.

Also, can you add the test_Hessian to one of the existing test scripts such that it will always run?

src/xSTIR/cSTIR/cstir.cpp

src/xSTIR/pSTIR/STIR.py

examples/Python/PET/osem_reconstruction.py

KrisThielemans · 2024-05-16T00:51:02Z

Also, please comment out the Coveralls lines
https://github.com/SyneRBI/SIRF/blob/d54c4550d1533132c6a281b6cf2159c8fd22f45c/.github/workflows/build-test.yml#L133C4-L147C32
You probably have to preserve indentation. Or leave that for a separate PR

KrisThielemans · 2024-05-16T11:38:59Z

@evgueni-ovtchinnikov test_Hessian isn't a great test, as it just compares the norms, while it should compare the difference.

I'm modifying it and putting it in test_ObjectiveFunction.py

  def test_Hessian(self, subset=-1, eps=1e-3):
        """Checks that grad(x + dx) - grad(x) is close to H(x)*dx
        """
        x = self.image
        dx = x.clone()
        dx *= eps/dx.norm()
        dx += eps/2
        y = x + dx
        gx = self.obj_fun.gradient(x, subset)
        gy = self.obj_fun.gradient(y, subset)
        dg = gy - gx
        Hdx = self.obj_fun.multiply_with_Hessian(x, dx, subset)
        norm = dg.norm()
        q = (dg - Hdx).norm()/dg.norm()
        print('norm of grad(x + dx) - grad(x): %f' % dg.norm())
        print('norm of H(x)*dx: %f' % Hdx.norm())
        print('relative difference: %f' % q)
        numpy.testing.assert_array_equal(q,0)

However, I had to add a constant term to the acq_model as otherwise the test fails (essentially because the Hessian is ill-defined otherwise). I'll commit this soon.

KrisThielemans · 2024-05-16T12:21:25Z

@evgueni-ovtchinnikov I'm now adding a similar test to the priors. (Had to delete a subset in the call). However, I'm confused by results. Does get_gradient(x) modify the result? Looks like it for priors.

evgueni-ovtchinnikov · 2024-05-16T13:07:35Z

@KrisThielemans

Does get_gradient(x) modify the result? Looks like it for priors.

what result?

KrisThielemans · 2024-05-16T14:11:01Z

got distracted... Sorry,

I did push the test for the LogLikelihood. Works for me.

KrisThielemans · 2024-05-16T14:37:15Z

@KrisThielemans

Does get_gradient(x) modify the result? Looks like it for priors.

what result?

ignore that. The loop in the test is a bit hard to understand... Update coming soon.

also remove test_Hessian from the objective function

KrisThielemans · 2024-05-16T14:45:37Z

ok. Should be all done now. Can you add a line to CHANGES.md and merge?

evgueni-ovtchinnikov added 5 commits May 2, 2024 14:21

implemented prior.accumulate_Hessian_times_input()

6816fc0

implemented return of Hessian*x via out= in Prior

956c498

merged master

2390735

implemented ObjectiveFunction.accumulate_Hessian_times_input (not tes…

99ab3c6

…ted)

implemented set_time_interval for listmode objective function (not te…

c14da39

…sted)

KrisThielemans requested changes May 10, 2024

View reviewed changes

KrisThielemans mentioned this pull request May 10, 2024

list-mode obj_fun.get_subset_sensitivity(0) returns None #1251

Closed

evgueni-ovtchinnikov added 2 commits May 10, 2024 16:46

used auto and corrected inheritance of listmode objective function in…

f5ff1ba

… Python

fixed typo in listmode objective function class in STIR.py

66c35c7

use references in Hessian

14bf4f6

Also use auto more

evgueni-ovtchinnikov added 2 commits May 15, 2024 11:12

implemented multiply_with_Hessian, checked by simple regression test

0bb6619

see osem_reconstruction.py for a simple regression test for the Hessian multiplication

attended to Codacy issue

136add7

gschramm mentioned this pull request May 15, 2024

listmode DL notebooks for PSMR 2024 SyneRBI/SIRF-Exercises#225

Merged

added objective function method that checks that grad(x+dx)-grad(x)~H…

7f24a44

…(x)dx

evgueni-ovtchinnikov marked this pull request as ready for review May 15, 2024 16:37

attended to Codacy issues

4807871

KrisThielemans requested changes May 16, 2024

View reviewed changes

resolved reviewer's issues

3b2dffe

updated Hessian multiplication for priors

608631c

evgueni-ovtchinnikov and others added 2 commits May 16, 2024 13:20

corrected Prior.multiply_with_Hessian()

16547a5

add test on Hessian to CI

32167f6

add test for Hessian of prior

2dbcacc

also remove test_Hessian from the objective function

evgueni-ovtchinnikov added 2 commits May 16, 2024 15:04

attended to Codacy issues

091de32

updated CHANGES.md

0531235

evgueni-ovtchinnikov merged commit 12949b4 into master May 16, 2024
6 checks passed

evgueni-ovtchinnikov deleted the acc-hess branch May 16, 2024 15:34

evgueni-ovtchinnikov mentioned this pull request Jun 27, 2024

time frame handling in listmode reconstruction #1249

Closed

Computing Hessian products with image data objects and related stuff #1253

Computing Hessian products with image data objects and related stuff #1253

Conversation

evgueni-ovtchinnikov commented May 10, 2024 • edited Loading

Changes in this pull request

Testing performed

Related issues

Checklist before requesting a review

Contribution Notes

KrisThielemans May 10, 2024

Choose a reason for hiding this comment

evgueni-ovtchinnikov May 10, 2024

Choose a reason for hiding this comment

evgueni-ovtchinnikov May 13, 2024

Choose a reason for hiding this comment

KrisThielemans May 15, 2024

Choose a reason for hiding this comment

KrisThielemans May 10, 2024

Choose a reason for hiding this comment

evgueni-ovtchinnikov May 10, 2024

Choose a reason for hiding this comment

KrisThielemans commented May 10, 2024

evgueni-ovtchinnikov commented May 10, 2024

gschramm commented May 12, 2024 • edited Loading

KrisThielemans commented May 12, 2024

gschramm commented May 12, 2024

KrisThielemans commented May 12, 2024

gschramm commented May 12, 2024

KrisThielemans commented May 12, 2024

gschramm commented May 13, 2024

KrisThielemans commented May 13, 2024

gschramm commented May 13, 2024

KrisThielemans commented May 13, 2024

gschramm commented May 13, 2024

KrisThielemans commented May 13, 2024

KrisThielemans commented May 13, 2024

KrisThielemans commented May 13, 2024

evgueni-ovtchinnikov commented May 13, 2024

gschramm commented May 13, 2024 • edited Loading

KrisThielemans commented May 13, 2024

KrisThielemans commented May 13, 2024

gschramm commented May 15, 2024

KrisThielemans commented May 15, 2024

KrisThielemans commented May 15, 2024

KrisThielemans commented May 15, 2024

evgueni-ovtchinnikov commented May 15, 2024

evgueni-ovtchinnikov commented May 15, 2024

KrisThielemans left a comment

Choose a reason for hiding this comment

KrisThielemans commented May 16, 2024

KrisThielemans commented May 16, 2024

KrisThielemans commented May 16, 2024

evgueni-ovtchinnikov commented May 16, 2024 • edited Loading

KrisThielemans commented May 16, 2024 • edited Loading

KrisThielemans commented May 16, 2024

KrisThielemans commented May 16, 2024

evgueni-ovtchinnikov commented May 10, 2024 •

edited

Loading

gschramm commented May 12, 2024 •

edited

Loading

gschramm commented May 13, 2024 •

edited

Loading

evgueni-ovtchinnikov commented May 16, 2024 •

edited

Loading

KrisThielemans commented May 16, 2024 •

edited

Loading