From 728efaf8e9ec316d09ed5c922d215ae5d22b6063 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca <14138589+paskino@users.noreply.github.com> Date: Thu, 22 Aug 2024 11:12:03 +0100 Subject: [PATCH 1/5] remove multiple exit from prange closes #1900 Signed-off-by: Edoardo Pasca <14138589+paskino@users.noreply.github.com> --- Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py index 9fa890a6f0..b990c0c638 100644 --- a/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py @@ -354,7 +354,7 @@ def kl_div(x, y, eta): accumulator[get_thread_id()] += Y else: # out.flat[i] = numpy.inf - return numpy.inf + accumulator[get_thread_id()] = numpy.inf return sum(accumulator) @njit(parallel=True) @@ -372,7 +372,7 @@ def kl_div_mask(x, y, eta, mask): accumulator[get_thread_id()] += Y else: # out.flat[i] = numpy.inf - return numpy.inf + accumulator[get_thread_id()] = numpy.inf return sum(accumulator) # convex conjugate From 07523228c3524662ac5b41c6aaa014bb77d3fcf2 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca <14138589+paskino@users.noreply.github.com> Date: Tue, 24 Sep 2024 15:02:49 +0100 Subject: [PATCH 2/5] Added description of kl_div fix Signed-off-by: Edoardo Pasca <14138589+paskino@users.noreply.github.com> --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index a87ce678fc..9118372719 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ - Use ravel instead of flat in KullbackLeibler numba backend (#1874) - Upgrade Python wrapper (#1873, #1875) - Add checks on out argument passed to processors to ensure corrrect dtype and size (#1805) + - Removed multiple exits from numba implementation of KullbackLeibler divergence (#1901) - Testing: - New unit tests for operators and functions to check for in place errors and the behaviour of `out` (#1805) - Bug fixes: From 5fe4179ffbb3d6f4519291e547cd6ac88e2aecb1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca <14138589+paskino@users.noreply.github.com> Date: Tue, 1 Oct 2024 08:38:30 +0100 Subject: [PATCH 3/5] skips computation of KL divergence if inf was already computed --- .../optimisation/functions/KullbackLeibler.py | 42 +++++++++++-------- 1 file changed, 24 insertions(+), 18 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py index b990c0c638..cc1f394cf7 100644 --- a/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py @@ -343,25 +343,9 @@ def kl_gradient_mask(x, b, out, eta, mask): @njit(parallel=True) def kl_div(x, y, eta): accumulator = numpy.zeros(get_num_threads(), dtype=numpy.float64) + has_inf = numpy.zeros(get_num_threads(), dtype=numpy.int8) for i in prange(x.size): - X = x.flat[i] - Y = y.flat[i] + eta.flat[i] - if X > 0 and Y > 0: - # out.flat[i] = X * numpy.log(X/Y) - X + Y - accumulator[get_thread_id()] += X * numpy.log(X/Y) - X + Y - elif X == 0 and Y >= 0: - # out.flat[i] = Y - accumulator[get_thread_id()] += Y - else: - # out.flat[i] = numpy.inf - accumulator[get_thread_id()] = numpy.inf - return sum(accumulator) - - @njit(parallel=True) - def kl_div_mask(x, y, eta, mask): - accumulator = numpy.zeros(get_num_threads(), dtype=numpy.float64) - for i in prange(x.size): - if mask.flat[i] > 0: + if has_inf[get_thread_id()] == 0: X = x.flat[i] Y = y.flat[i] + eta.flat[i] if X > 0 and Y > 0: @@ -373,6 +357,28 @@ def kl_div_mask(x, y, eta, mask): else: # out.flat[i] = numpy.inf accumulator[get_thread_id()] = numpy.inf + has_inf[get_thread_id()] = 1 + return sum(accumulator) + + @njit(parallel=True) + def kl_div_mask(x, y, eta, mask): + accumulator = numpy.zeros(get_num_threads(), dtype=numpy.float64) + has_inf = numpy.zeros(get_num_threads(), dtype=numpy.int8) + for i in prange(x.size): + if has_inf[get_thread_id()] == 0: + if mask.flat[i] > 0: + X = x.flat[i] + Y = y.flat[i] + eta.flat[i] + if X > 0 and Y > 0: + # out.flat[i] = X * numpy.log(X/Y) - X + Y + accumulator[get_thread_id()] += X * numpy.log(X/Y) - X + Y + elif X == 0 and Y >= 0: + # out.flat[i] = Y + accumulator[get_thread_id()] += Y + else: + # out.flat[i] = numpy.inf + accumulator[get_thread_id()] = numpy.inf + has_inf[get_thread_id()] = 1 return sum(accumulator) # convex conjugate From 9c23739ab01c5811b9fa4f4abb67177cbe692c77 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca <14138589+paskino@users.noreply.github.com> Date: Wed, 2 Oct 2024 09:26:42 +0100 Subject: [PATCH 4/5] use single int for stopping threads and add unit test --- .../optimisation/functions/KullbackLeibler.py | 12 ++++++------ .../Python/test/test_function_KullbackLeibler.py | 16 ++++++++++++++++ 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py index cc1f394cf7..3a84fcc8d9 100644 --- a/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py +++ b/Wrappers/Python/cil/optimisation/functions/KullbackLeibler.py @@ -343,9 +343,9 @@ def kl_gradient_mask(x, b, out, eta, mask): @njit(parallel=True) def kl_div(x, y, eta): accumulator = numpy.zeros(get_num_threads(), dtype=numpy.float64) - has_inf = numpy.zeros(get_num_threads(), dtype=numpy.int8) + has_inf = 0 for i in prange(x.size): - if has_inf[get_thread_id()] == 0: + if has_inf == 0: X = x.flat[i] Y = y.flat[i] + eta.flat[i] if X > 0 and Y > 0: @@ -357,15 +357,15 @@ def kl_div(x, y, eta): else: # out.flat[i] = numpy.inf accumulator[get_thread_id()] = numpy.inf - has_inf[get_thread_id()] = 1 + has_inf = 1 return sum(accumulator) @njit(parallel=True) def kl_div_mask(x, y, eta, mask): accumulator = numpy.zeros(get_num_threads(), dtype=numpy.float64) - has_inf = numpy.zeros(get_num_threads(), dtype=numpy.int8) + has_inf = 0 for i in prange(x.size): - if has_inf[get_thread_id()] == 0: + if has_inf == 0: if mask.flat[i] > 0: X = x.flat[i] Y = y.flat[i] + eta.flat[i] @@ -378,7 +378,7 @@ def kl_div_mask(x, y, eta, mask): else: # out.flat[i] = numpy.inf accumulator[get_thread_id()] = numpy.inf - has_inf[get_thread_id()] = 1 + has_inf = 1 return sum(accumulator) # convex conjugate diff --git a/Wrappers/Python/test/test_function_KullbackLeibler.py b/Wrappers/Python/test/test_function_KullbackLeibler.py index e5e9dc9ed4..a71ccfea5c 100644 --- a/Wrappers/Python/test/test_function_KullbackLeibler.py +++ b/Wrappers/Python/test/test_function_KullbackLeibler.py @@ -173,6 +173,22 @@ def test_KullbackLeibler_numba_call(self): numpy.testing.assert_allclose(f(u1), f_np(u1), rtol=1e-5) + def test_KullbackLeibler_numba_kl_div_has_inf(self): + x = self.u1 * self.mask + x *= -1 + + from cil.optimisation.functions.KullbackLeibler import kl_div + + numpy.testing.assert_equal(kl_div(x.array, self.b1.array, self.eta.array), numpy.inf) + + def test_KullbackLeibler_numba_kl_div_mask_has_inf(self): + x = self.u1 * self.mask + x *= -1 + + from cil.optimisation.functions.KullbackLeibler import kl_div_mask + + numpy.testing.assert_equal(kl_div_mask(x.array, self.b1.array, self.eta.array, self.mask.array), numpy.inf) + def test_KullbackLeibler_numba_call_mask(self): f = self.f f_np = self.f_np From f16cce6fb1e74b80b3386afb2ae1b6d759779c26 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca <14138589+paskino@users.noreply.github.com> Date: Wed, 16 Oct 2024 13:02:45 +0100 Subject: [PATCH 5/5] updated changelog --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 88f53e0df2..39894b7278 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,8 @@ * 24.x.x - Bug fixes: - Fix bug with 'median' and 'mean' methods in Masker averaging over the wrong axes. +- Enhancements: + - Removed multiple exits from numba implementation of KullbackLeibler divergence (#1901) * 24.2.0 - New Features: @@ -22,7 +24,6 @@ - Allow the SIRT algorithm to take `initial=None` (#1906) - Add checks on equality method of `AcquisitionData` and `ImageData` for equality of data type and geometry (#1919) - Add check on equality method of `AcquisitionGeometry` for equality of dimension labels (#1919) - - Removed multiple exits from numba implementation of KullbackLeibler divergence (#1901) - Testing: - New unit tests for operators and functions to check for in place errors and the behaviour of `out` (#1805) - Updates in SPDHG vs PDHG unit test to reduce test time and adjustments to parameters (#1898)