xapyb pet #895

gfardell · 2021-04-20T11:13:37Z

There is no xapyb where a and b are mixed scalar and vector, therefore I've implemented this in python with the vector multiplication first before the scalar version of axpby is called. If we want to implement this at a lower level I am happy to, but it'll need it to be added to STIR too.
I have not implemented anything in matlab, I don't have a matlab to test the code and don't know matlab well at all, so it would be good if someone else could handle this
I've added unit tests at the python level
I haven't implemented the vector xapyb versions for gadgetron data containers
Where python tries to call the vector version and this is not implemented the error is caught and python is used for the calculation instead however it still prints:

File: /home/sirfuser/devel/install/python/sirf/SIRF.py
Line: 254
check_status found the following message sent from the engine:

I'm not sure if there's a way to handle this try/except better.

KrisThielemans

this was quite a lot of work! thanks

KrisThielemans · 2021-04-20T13:54:12Z

src/Registration/cReg/NiftiImageData.cpp

+    const NiftiImageData<dataType>& a = dynamic_cast<const NiftiImageData<dataType>&>(a_a);
+    const NiftiImageData<dataType>& b = dynamic_cast<const NiftiImageData<dataType>&>(a_b);
+    const NiftiImageData<dataType>& x = dynamic_cast<const NiftiImageData<dataType>&>(a_x);
+    const NiftiImageData<dataType>& y = dynamic_cast<const NiftiImageData<dataType>&>(a_y);


could potentially throw unexpectedly if the type isn't ok. This is fine, but will confuse the user (certainly so a Python user).

should we cope with cases where the data-type isn't NiftiImageData? Might be quite hard. But if we don't, I think we should put the whole thing in a try catch block with some sensible error.

Also, in future, this can be shorter as

auto& a = dynamic_cast<const NiftiImageData<dataType>&>(a_a);

At the python end sapyb for vectors is in a try/except so it will default to python if there's an issue.

It does seem a bit strange restricting a and b to the same types, but I don't know how data is organised so whether a numpy array (for example) is compatible with nifty or stir arrays.

ok. Nevertheless, it's the C++ library that should give sensible errors as much as we can, not the "client" Python/MATLAB

src/Registration/cReg/NiftiImageData.cpp

KrisThielemans · 2021-04-20T13:56:09Z

src/common/SIRF.py

@@ -202,6 +202,8 @@ def add(self, other, out=None):
        return z
    def axpby(self, a, b, y, out=None, **kwargs):
        '''
+        This function has been deprecated and will be removed in a future release. Please use `sapyb` instead.


better to use https://python-deprecated.readthedocs.io/en/latest/tutorial.html?

I'm happy to use this. But I'll need to add the Python-Deprecated to the requirements, could you point out where these live?

should be there already
https://github.com/SyneRBI/SIRF/blob/master/requirements.txt#L7

evgueni-ovtchinnikov · 2021-04-20T14:06:00Z

@KrisThielemans

#define THROW(msg) throw LocalisedException(msg, __FILE__, __LINE__)
#define ASSERT(condition, msg) if (!(condition)) THROW(msg)

KrisThielemans · 2021-04-20T14:10:37Z

ok!

gfardell · 2021-04-20T15:24:40Z

this was quite a lot of work! thanks

400+ lines of code does seem a little excessive for basically some wrappers, but at least I understand SIRF a little now!

src/common/SIRF.py

evgueni-ovtchinnikov · 2021-04-20T14:45:58Z

src/common/SIRF.py

+                    tmp = self.multiply(a)
+                    y.multiply(b, out=z)
+                    z.add(tmp, out=z)


should better check the error message, and if it is "xapyb has not been implemented yet" then return NotImplemented, otherwise throw (something else happened)

Frankly, I feel uneasy about your trying to cover lots of cases (axpby, xapby and essentially also axpyb and xapby + self variants) by one method - are all of them really needed? If they all are, I would rather have had several methods.

besides, we normally try to do as little as possible in Python and Matlab interfaces, if only to reduce duplication

axpby previously covered vector and scalar methods, but with python computing the vector case, so I don't think it's covering more cases. But yes, it is a busy method.

In terms of the users passing a data-container to out that is also one of the inputs, I guess we'd have to check on the ID to stop this, but as the operation is pixel-wise anyway I assumed it would be safe to allow. It doesn't change the implementation apart from the python code order avoids an intermediary step writing to out.

I think CIL currently wraps the C++ code with ctypes and passes the type of a and b as parameters. I guess this would be an option to simplify the python and matlab ends. but sapyb needs to be consistent in behaviour across SIRF and CIL, I'm not sure we'd want to split it in to sapyb_vv or sapyb_ss at the python end as there's algorithms where a and b could be either.

using input container as out is ok - see e.g. __iadd__.

I still believe unnecessary copies can be avoided, but let us leave this for now - after all, as I mentioned, in SIRF we try to do all computation in C++, so sapyb will slim down eventually.

The only unnecessary copy is when a and b are of different types, and the user has passed out, as it will pre-calculate the vector * vector side. Is this where you mean?

There is a specific case if the users passes y and out as the same object, where a is a vector and b a scalar. Then an intermediate step x.multiply(a, out=out) would overwrite y before it was used. So without checking on the datacontainer IDs it seems safer that where it has to calculate in steps it writes initially for a temporary container rather than out. Of course if it was implemented in the C++ code it would be pixel-wise and not an issue.

KrisThielemans · 2021-04-20T16:04:04Z

Turns out we have a DEPRECATED macro already in

SIRF/src/iUtilities/include/sirf/iUtilities/iutilities.h

Line 24 in 843eb25

// Deprecation function. With C++14, could use [[deprecated("some message")]]

I like it better that it has its own file though. I think SIRF_DEPRECATED might be better than DEPRECATED, but that seems to have to be for a different PR

KrisThielemans · 2021-04-21T12:17:30Z

@gfardell @evgueni-ovtchinnikov I have no clear feeling if this is feasible to merge in the next few days, or if we better postpone it (which would be fine). can you comment?

gfardell · 2021-04-21T12:36:12Z

@gfardell @evgueni-ovtchinnikov I have no clear feeling if this is feasible to merge in the next few days, or if we better postpone it (which would be fine). can you comment?

I am updating it now and don't imagine it will take too long for it's current functionality, however if we want to delay it until the gadegtron xapyb is implemented I guess it would be more complete.

evgueni-ovtchinnikov · 2021-04-21T12:48:34Z

@KrisThielemans I would postpone it - for one thing, I believe we should move the computation in Python sapyb to C++ and do Matlab and MR before merge. (Now that Johannes took over Resampler from me, I can take over xapyb from Gemma :-))

gfardell · 2021-04-22T12:33:57Z

I can't figure out why this stopped building. I can build and run the tests locally still.

KrisThielemans · 2021-04-22T13:10:58Z

As the error is with recent changes in master, I think it is best to merge master on here again as the DEPRECATED stuff changed there. Presumably might give a conflict (although github doesn't think so).

gfardell · 2021-05-04T09:24:05Z

I've made a stab at the matlab wrappers but I haven't managed to get a license from SCD yet to run and test it. I've also not used matlab much so have no clue on error handling, when I get a license I'll be able to figure that out.

evgueni-ovtchinnikov · 2021-05-04T13:55:57Z

src/common/+sirf/+SIRF/DataContainer.m

+                tmp = x.same_object();
+
+                if a_scalar
+                    tmp = b.dot(y);


should be b.times(y), dot computes scalar product

evgueni-ovtchinnikov · 2021-05-04T13:56:12Z

src/common/+sirf/+SIRF/DataContainer.m

+                    z.handle_ = calllib('msirf', 'mSIRF_axpby', ...
+                        ptr_a, x.handle_, 1.0, tmp.handle_);
+                else
+                    tmp = a.dot(x);


evgueni-ovtchinnikov · 2021-05-04T14:05:10Z

src/common/+sirf/+SIRF/DataContainer.m

+                z.handle_ = calllib('msirf', 'mSIRF_axpby', ...
+                    ptr_a, x.handle_, ptr_b, y.handle_);


if neither a nor b is scalar, should call mSIRF_xapyb_vv

evgueni-ovtchinnikov · 2021-05-05T09:00:29Z

src/common/csirf.cpp

+void*
+cSIRF_xapyb_ss(
+const void* ptr_x, const void* ptr_a,
+const void* ptr_y, const void* ptr_b
+) {
+	try {
+		DataContainer& x =
+			objectFromHandle<DataContainer >(ptr_x);
+		DataContainer& y =
+			objectFromHandle<DataContainer >(ptr_y);
+		void* h = x.new_data_container_handle();
+		DataContainer& z = objectFromHandle<DataContainer>(h);
+		z.xapyb(x, ptr_a, y, ptr_b);
+		return h;
+	}
+	CATCH;
+}
+
+extern "C"
+void*
+cSIRF_xapyb_vv(
+const void* ptr_x, const void* ptr_a,
+const void* ptr_y, const void* ptr_b
+) {
+	try {
+		DataContainer& x =
+			objectFromHandle<DataContainer >(ptr_x);
+		DataContainer& y =


these two functions have substantial overlap - why not have one cSIRF_xapyb with extra const char* args parameter that may have values "ss" or "vv" (and potentially other values, making it future proof)?

this would also solve lesser problem of breaking our naming scheme cSIRF_nameOfFunctionInCamelCase

evgueni-ovtchinnikov · 2021-05-05T09:05:51Z

src/common/include/sirf/common/csirf.h

+void* cSIRF_xapyb_ss(const void* ptr_x, const PTR_FLOAT ptr_a,
+	const void* ptr_y, const PTR_FLOAT ptr_b);
+void* cSIRF_xapyb_ss_Alt(const void* ptr_x, const PTR_FLOAT ptr_a,
+	const void* ptr_y, const PTR_FLOAT ptr_b, void* ptr_z);
+void* cSIRF_xapyb_vv(const void* ptr_x, const void*  ptr_a,
+	const void* ptr_y, const void* ptr_b);
+void* cSIRF_xapyb_vv_Alt(const void* ptr_x, const void*  ptr_a,
+	const void* ptr_y, const void* ptr_b, void* ptr_z);	


see the above suggestion to replace these four functions with two functions with an extra const char* args parameter

KrisThielemans · 2021-05-20T09:21:46Z

merged via another branch

gfardell added 7 commits April 14, 2021 16:00

add xapyb to sirf datacontainers

5f7839c

c wrapeprs and python implementation

3ccb474

throw where not implemented

0a299de

add deprecation notice

7abca0e

bug fixes

546966c

clean up

98da8bb

bugfix

6c15dce

gfardell requested review from evgueni-ovtchinnikov, paskino and KrisThielemans and removed request for evgueni-ovtchinnikov April 20, 2021 11:13

added nify image xapyb for vectors

e39cb70

KrisThielemans linked an issue Apr 20, 2021 that may be closed by this pull request

axpby renaming to xapby #834

Closed

KrisThielemans requested changes Apr 20, 2021

View reviewed changes

evgueni-ovtchinnikov requested changes Apr 20, 2021

View reviewed changes

KrisThielemans mentioned this pull request Apr 21, 2021

rename DEPRECATED to SIRF_DEPRECATED #900

Closed

gfardell added 3 commits April 21, 2021 12:55

Updated xapby methods called

dbb6e6c

clean up

11a6159

deprecate axpby

bb122b6

KrisThielemans added this to the v3.1 milestone Apr 21, 2021

gfardell added 2 commits April 21, 2021 16:12

nifty try-catch

9a1d473

white space fix [CI SKIP]

d93f1da

gfardell added 4 commits April 22, 2021 13:58

Merge remote-tracking branch 'origin/master' into xapyb_pet

d7e02e2

remove duplicated deprecate macro

4d7fa14

remove redundat deprecate macro

0055848

matlab wrapper initial

eb3bde0

evgueni-ovtchinnikov requested changes May 5, 2021

View reviewed changes

KrisThielemans closed this May 20, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

xapyb pet #895

xapyb pet #895

gfardell commented Apr 20, 2021 •

edited

Loading

KrisThielemans left a comment

KrisThielemans Apr 20, 2021

gfardell Apr 20, 2021

KrisThielemans Apr 21, 2021

KrisThielemans Apr 20, 2021

gfardell Apr 21, 2021

KrisThielemans Apr 21, 2021

evgueni-ovtchinnikov commented Apr 20, 2021

KrisThielemans commented Apr 20, 2021

gfardell commented Apr 20, 2021

evgueni-ovtchinnikov Apr 20, 2021

evgueni-ovtchinnikov Apr 20, 2021

evgueni-ovtchinnikov Apr 20, 2021

gfardell Apr 20, 2021

evgueni-ovtchinnikov Apr 20, 2021

gfardell Apr 21, 2021

KrisThielemans commented Apr 20, 2021

KrisThielemans commented Apr 21, 2021

gfardell commented Apr 21, 2021

evgueni-ovtchinnikov commented Apr 21, 2021

gfardell commented Apr 22, 2021

KrisThielemans commented Apr 22, 2021

gfardell commented May 4, 2021

evgueni-ovtchinnikov May 4, 2021

evgueni-ovtchinnikov May 4, 2021

evgueni-ovtchinnikov May 4, 2021

evgueni-ovtchinnikov May 5, 2021

evgueni-ovtchinnikov May 5, 2021

KrisThielemans commented May 20, 2021

		z.handle_ = calllib('msirf', 'mSIRF_axpby', ...
		ptr_a, x.handle_, ptr_b, y.handle_);

xapyb pet #895

xapyb pet #895

Conversation

gfardell commented Apr 20, 2021 • edited Loading

KrisThielemans left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

evgueni-ovtchinnikov commented Apr 20, 2021

KrisThielemans commented Apr 20, 2021

gfardell commented Apr 20, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KrisThielemans commented Apr 20, 2021

KrisThielemans commented Apr 21, 2021

gfardell commented Apr 21, 2021

evgueni-ovtchinnikov commented Apr 21, 2021

gfardell commented Apr 22, 2021

KrisThielemans commented Apr 22, 2021

gfardell commented May 4, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KrisThielemans commented May 20, 2021

gfardell commented Apr 20, 2021 •

edited

Loading