summaryrefslogtreecommitdiff
path: root/silx/opencl/test
diff options
context:
space:
mode:
Diffstat (limited to 'silx/opencl/test')
-rw-r--r--silx/opencl/test/__init__.py7
-rw-r--r--silx/opencl/test/test_backprojection.py78
-rw-r--r--silx/opencl/test/test_kahan.py269
-rw-r--r--silx/opencl/test/test_stats.py114
4 files changed, 451 insertions, 17 deletions
diff --git a/silx/opencl/test/__init__.py b/silx/opencl/test/__init__.py
index 2fe88ea..f3121d5 100644
--- a/silx/opencl/test/__init__.py
+++ b/silx/opencl/test/__init__.py
@@ -24,7 +24,7 @@
__authors__ = ["J. Kieffer"]
__license__ = "MIT"
-__date__ = "17/10/2017"
+__date__ = "11/01/2019"
import os
import unittest
@@ -36,6 +36,9 @@ from . import test_linalg
from . import test_array_utils
from ..codec import test as test_codec
from . import test_image
+from . import test_kahan
+from . import test_stats
+
def suite():
test_suite = unittest.TestSuite()
@@ -47,6 +50,8 @@ def suite():
test_suite.addTests(test_array_utils.suite())
test_suite.addTests(test_codec.suite())
test_suite.addTests(test_image.suite())
+ test_suite.addTests(test_kahan.suite())
+ test_suite.addTests(test_stats.suite())
# Allow to remove sift from the project
test_base_dir = os.path.dirname(__file__)
sift_dir = os.path.join(test_base_dir, "..", "sift")
diff --git a/silx/opencl/test/test_backprojection.py b/silx/opencl/test/test_backprojection.py
index 70ce2ae..7ab416d 100644
--- a/silx/opencl/test/test_backprojection.py
+++ b/silx/opencl/test/test_backprojection.py
@@ -35,8 +35,9 @@ __date__ = "19/01/2018"
import time
import logging
-import numpy
+import numpy as np
import unittest
+from math import pi
try:
import mako
except ImportError:
@@ -44,6 +45,7 @@ except ImportError:
from ..common import ocl
if ocl:
from .. import backprojection
+ from ..sinofilter import compute_fourier_filter
from silx.test.utils import utilstest
logger = logging.getLogger(__name__)
@@ -55,7 +57,7 @@ def generate_coords(img_shp, center=None):
The zero is at the center of the image.
"""
l_r, l_c = float(img_shp[0]), float(img_shp[1])
- R, C = numpy.mgrid[:l_r, :l_c]
+ R, C = np.mgrid[:l_r, :l_c]
if center is None:
center0, center1 = l_r / 2., l_c / 2.
else:
@@ -71,7 +73,7 @@ def clip_circle(img, center=None, radius=None):
"""
R, C = generate_coords(img.shape, center)
M = R * R + C * C
- res = numpy.zeros_like(img)
+ res = np.zeros_like(img)
if radius is None:
radius = img.shape[0] / 2. - 1
mask = M < radius * radius
@@ -85,23 +87,29 @@ class TestFBP(unittest.TestCase):
def setUp(self):
if ocl is None:
return
- # ~ if sys.platform.startswith('darwin'):
- # ~ self.skipTest("Backprojection is not implemented on CPU for OS X yet")
self.getfiles()
self.fbp = backprojection.Backprojection(self.sino.shape, profile=True)
if self.fbp.compiletime_workgroup_size < 16 * 16:
- self.skipTest("Current implementation of OpenCL backprojection is not supported on this platform yet")
+ self.skipTest("Current implementation of OpenCL backprojection is "
+ "not supported on this platform yet")
+ # Astra does not use the same backprojector implementation.
+ # Therefore, we cannot expect results to be the "same" (up to float32
+ # numerical error)
+ self.tol = 5e-2
+ if self.fbp.is_cpu:
+ # Precision is less when using CPU
+ self.tol *= 2
def tearDown(self):
self.sino = None
-# self.fbp.log_profile()
+ # self.fbp.log_profile()
self.fbp = None
def getfiles(self):
# load sinogram of 512x512 MRI phantom
- self.sino = numpy.load(utilstest.getfile("sino500.npz"))["data"]
+ self.sino = np.load(utilstest.getfile("sino500.npz"))["data"]
# load reconstruction made with ASTRA FBP (with filter designed in spatial domain)
- self.reference_rec = numpy.load(utilstest.getfile("rec_astra_500.npz"))["data"]
+ self.reference_rec = np.load(utilstest.getfile("rec_astra_500.npz"))["data"]
def measure(self):
"Common measurement of timings"
@@ -124,8 +132,8 @@ class TestFBP(unittest.TestCase):
ref_clipped = clip_circle(self.reference_rec)
delta = abs(res_clipped - ref_clipped)
bad = delta > 1
-# numpy.save("/tmp/bad.npy", bad.astype(int))
- logger.debug("Absolute difference: %s with %s outlier pixels out of %s", delta.max(), bad.sum(), numpy.prod(bad.shape))
+ logger.debug("Absolute difference: %s with %s outlier pixels out of %s"
+ "", delta.max(), bad.sum(), np.prod(bad.shape))
return delta.max()
@unittest.skipUnless(ocl and mako, "pyopencl is missing")
@@ -143,21 +151,59 @@ class TestFBP(unittest.TestCase):
err = self.compare(res)
msg = str("Max error = %e" % err)
logger.info(msg)
- # TODO: cannot do better than 1e0 ?
- # The plain backprojection was much better, so it must be an issue in the filtering process
- self.assertTrue(err < 1., "Max error is too high")
+ self.assertTrue(err < self.tol, "Max error is too high")
+
# Test multiple reconstructions
# -----------------------------
- res0 = numpy.copy(res)
+ res0 = np.copy(res)
for i in range(10):
res = self.fbp.filtered_backprojection(self.sino)
- errmax = numpy.max(numpy.abs(res - res0))
+ errmax = np.max(np.abs(res - res0))
self.assertTrue(errmax < 1.e-6, "Max error is too high")
+ @unittest.skipUnless(ocl and mako, "pyopencl is missing")
+ def test_fbp_filters(self):
+ """
+ Test the different available filters of silx FBP.
+ """
+ avail_filters = [
+ "ramlak", "shepp-logan", "cosine", "hamming",
+ "hann"
+ ]
+ # Create a Dirac delta function at a single angle view.
+ # As the filters are radially invarant:
+ # - backprojection yields an image where each line is a Dirac.
+ # - FBP yields an image where each line is the spatial filter
+ # One can simply filter "dirac" without backprojecting it, but this
+ # test will also ensure that backprojection behaves well.
+ dirac = np.zeros_like(self.sino)
+ na, dw = dirac.shape
+ dirac[0, dw//2] = na / pi * 2
+
+ for filter_name in avail_filters:
+ B = backprojection.Backprojection(dirac.shape, filter_name=filter_name)
+ r = B(dirac)
+ # Check that radial invariance is kept
+ std0 = np.max(np.abs(np.std(r, axis=0)))
+ self.assertTrue(
+ std0 < 5.e-6,
+ "Something wrong with FBP(filter=%s)" % filter_name
+ )
+ # Check that the filter is retrieved
+ r_f = np.fft.fft(np.fft.fftshift(r[0])).real / 2. # filter factor
+ ref_filter_f = compute_fourier_filter(dw, filter_name)
+ errmax = np.max(np.abs(r_f - ref_filter_f))
+ logger.info("FBP filter %s: max error=%e" % (filter_name, errmax))
+ self.assertTrue(
+ errmax < 1.e-3,
+ "Something wrong with FBP(filter=%s)" % filter_name
+ )
+
def suite():
testSuite = unittest.TestSuite()
testSuite.addTest(TestFBP("test_fbp"))
+ testSuite.addTest(TestFBP("test_fbp_filters"))
return testSuite
diff --git a/silx/opencl/test/test_kahan.py b/silx/opencl/test/test_kahan.py
new file mode 100644
index 0000000..bb9ea3f
--- /dev/null
+++ b/silx/opencl/test/test_kahan.py
@@ -0,0 +1,269 @@
+#!/usr/bin/env python
+# coding: utf-8
+#
+# Project: Azimuthal integration
+# https://github.com/silx-kit/pyFAI
+#
+# Copyright (C) 2015-2019 European Synchrotron Radiation Facility, Grenoble, France
+#
+# Principal author: Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+
+"test suite for OpenCL code"
+
+from __future__ import absolute_import, division, print_function
+
+__author__ = "Jérôme Kieffer"
+__contact__ = "Jerome.Kieffer@ESRF.eu"
+__license__ = "MIT"
+__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
+__date__ = "30/01/2019"
+
+
+import unittest
+import numpy
+import logging
+import platform
+
+logger = logging.getLogger(__name__)
+try:
+ import pyopencl
+except ImportError as error:
+ logger.warning("OpenCL module (pyopencl) is not present, skip tests. %s.", error)
+ pyopencl = None
+
+from .. import ocl
+if ocl is not None:
+ from ..utils import read_cl_file
+ from .. import pyopencl
+ import pyopencl.array
+from ...test.utils import test_options
+
+
+class TestKahan(unittest.TestCase):
+ """
+ Test the kernels for compensated math in OpenCL
+ """
+
+ @classmethod
+ def setUpClass(cls):
+ if not test_options.WITH_OPENCL_TEST:
+ raise unittest.SkipTest("User request to skip OpenCL tests")
+ if pyopencl is None or ocl is None:
+ raise unittest.SkipTest("OpenCL module (pyopencl) is not present or no device available")
+
+ cls.ctx = ocl.create_context(devicetype="GPU")
+ cls.queue = pyopencl.CommandQueue(cls.ctx, properties=pyopencl.command_queue_properties.PROFILING_ENABLE)
+
+ # this is running 32 bits OpenCL woth POCL
+ if (platform.machine() in ("i386", "i686", "x86_64") and (tuple.__itemsize__ == 4) and
+ cls.ctx.devices[0].platform.name == 'Portable Computing Language'):
+ cls.args = "-DX87_VOLATILE=volatile"
+ else:
+ cls.args = ""
+
+ @classmethod
+ def tearDownClass(cls):
+ cls.queue = None
+ cls.ctx = None
+
+ @staticmethod
+ def dummy_sum(ary, dtype=None):
+ "perform the actual sum in a dummy way "
+ if dtype is None:
+ dtype = ary.dtype.type
+ sum_ = dtype(0)
+ for i in ary:
+ sum_ += i
+ return sum_
+
+ def test_kahan(self):
+ # simple test
+ N = 26
+ data = (1 << (N - 1 - numpy.arange(N))).astype(numpy.float32)
+
+ ref64 = numpy.sum(data, dtype=numpy.float64)
+ ref32 = self.dummy_sum(data)
+ if (ref64 == ref32):
+ logger.warning("Kahan: invalid tests as float32 provides the same result as float64")
+ # Dummy kernel to evaluate
+ src = """
+ kernel void summation(global float* data,
+ int size,
+ global float* result)
+ {
+ float2 acc = (float2)(0.0f, 0.0f);
+ for (int i=0; i<size; i++)
+ {
+ acc = kahan_sum(acc, data[i]);
+ }
+ result[0] = acc.s0;
+ result[1] = acc.s1;
+ }
+ """
+ prg = pyopencl.Program(self.ctx, read_cl_file("kahan.cl") + src).build(self.args)
+ ones_d = pyopencl.array.to_device(self.queue, data)
+ res_d = pyopencl.array.zeros(self.queue, 2, numpy.float32)
+ evt = prg.summation(self.queue, (1,), (1,), ones_d.data, numpy.int32(N), res_d.data)
+ evt.wait()
+ res = res_d.get().sum(dtype=numpy.float64)
+ self.assertEqual(ref64, res, "test_kahan")
+
+ def test_dot16(self):
+ # simple test
+ N = 16
+ data = (1 << (N - 1 - numpy.arange(N))).astype(numpy.float32)
+
+ ref64 = numpy.dot(data.astype(numpy.float64), data.astype(numpy.float64))
+ ref32 = numpy.dot(data, data)
+ if (ref64 == ref32):
+ logger.warning("dot16: invalid tests as float32 provides the same result as float64")
+ # Dummy kernel to evaluate
+ src = """
+ kernel void test_dot16(global float* data,
+ int size,
+ global float* result)
+ {
+ float2 acc = (float2)(0.0f, 0.0f);
+ float16 data16 = (float16) (data[0],data[1],data[2],data[3],data[4],
+ data[5],data[6],data[7],data[8],data[9],
+ data[10],data[11],data[12],data[13],data[14],data[15]);
+ acc = comp_dot16(data16, data16);
+ result[0] = acc.s0;
+ result[1] = acc.s1;
+ }
+
+ kernel void test_dot8(global float* data,
+ int size,
+ global float* result)
+ {
+ float2 acc = (float2)(0.0f, 0.0f);
+ float8 data0 = (float8) (data[0],data[2],data[4],data[6],data[8],data[10],data[12],data[14]);
+ float8 data1 = (float8) (data[1],data[3],data[5],data[7],data[9],data[11],data[13],data[15]);
+ acc = comp_dot8(data0, data1);
+ result[0] = acc.s0;
+ result[1] = acc.s1;
+ }
+
+ kernel void test_dot4(global float* data,
+ int size,
+ global float* result)
+ {
+ float2 acc = (float2)(0.0f, 0.0f);
+ float4 data0 = (float4) (data[0],data[4],data[8],data[12]);
+ float4 data1 = (float4) (data[3],data[7],data[11],data[15]);
+ acc = comp_dot4(data0, data1);
+ result[0] = acc.s0;
+ result[1] = acc.s1;
+ }
+
+ kernel void test_dot3(global float* data,
+ int size,
+ global float* result)
+ {
+ float2 acc = (float2)(0.0f, 0.0f);
+ float3 data0 = (float3) (data[0],data[4],data[12]);
+ float3 data1 = (float3) (data[3],data[11],data[15]);
+ acc = comp_dot3(data0, data1);
+ result[0] = acc.s0;
+ result[1] = acc.s1;
+ }
+
+ kernel void test_dot2(global float* data,
+ int size,
+ global float* result)
+ {
+ float2 acc = (float2)(0.0f, 0.0f);
+ float2 data0 = (float2) (data[0],data[14]);
+ float2 data1 = (float2) (data[1],data[15]);
+ acc = comp_dot2(data0, data1);
+ result[0] = acc.s0;
+ result[1] = acc.s1;
+ }
+
+ """
+
+ prg = pyopencl.Program(self.ctx, read_cl_file("kahan.cl") + src).build(self.args)
+ ones_d = pyopencl.array.to_device(self.queue, data)
+ res_d = pyopencl.array.zeros(self.queue, 2, numpy.float32)
+ evt = prg.test_dot16(self.queue, (1,), (1,), ones_d.data, numpy.int32(N), res_d.data)
+ evt.wait()
+ res = res_d.get().sum(dtype="float64")
+ self.assertEqual(ref64, res, "test_dot16")
+
+ res_d.fill(0)
+ data0 = data[0::2]
+ data1 = data[1::2]
+ ref64 = numpy.dot(data0.astype(numpy.float64), data1.astype(numpy.float64))
+ ref32 = numpy.dot(data0, data1)
+ if (ref64 == ref32):
+ logger.warning("dot8: invalid tests as float32 provides the same result as float64")
+ evt = prg.test_dot8(self.queue, (1,), (1,), ones_d.data, numpy.int32(N), res_d.data)
+ evt.wait()
+ res = res_d.get().sum(dtype="float64")
+ self.assertEqual(ref64, res, "test_dot8")
+
+ res_d.fill(0)
+ data0 = data[0::4]
+ data1 = data[3::4]
+ ref64 = numpy.dot(data0.astype(numpy.float64), data1.astype(numpy.float64))
+ ref32 = numpy.dot(data0, data1)
+ if (ref64 == ref32):
+ logger.warning("dot4: invalid tests as float32 provides the same result as float64")
+ evt = prg.test_dot4(self.queue, (1,), (1,), ones_d.data, numpy.int32(N), res_d.data)
+ evt.wait()
+ res = res_d.get().sum(dtype="float64")
+ self.assertEqual(ref64, res, "test_dot4")
+
+ res_d.fill(0)
+ data0 = numpy.array([data[0], data[4], data[12]])
+ data1 = numpy.array([data[3], data[11], data[15]])
+ ref64 = numpy.dot(data0.astype(numpy.float64), data1.astype(numpy.float64))
+ ref32 = numpy.dot(data0, data1)
+ if (ref64 == ref32):
+ logger.warning("dot3: invalid tests as float32 provides the same result as float64")
+ evt = prg.test_dot3(self.queue, (1,), (1,), ones_d.data, numpy.int32(N), res_d.data)
+ evt.wait()
+ res = res_d.get().sum(dtype="float64")
+ self.assertEqual(ref64, res, "test_dot3")
+
+ res_d.fill(0)
+ data0 = numpy.array([data[0], data[14]])
+ data1 = numpy.array([data[1], data[15]])
+ ref64 = numpy.dot(data0.astype(numpy.float64), data1.astype(numpy.float64))
+ ref32 = numpy.dot(data0, data1)
+ if (ref64 == ref32):
+ logger.warning("dot2: invalid tests as float32 provides the same result as float64")
+ evt = prg.test_dot2(self.queue, (1,), (1,), ones_d.data, numpy.int32(N), res_d.data)
+ evt.wait()
+ res = res_d.get().sum(dtype="float64")
+ self.assertEqual(ref64, res, "test_dot2")
+
+
+def suite():
+ testsuite = unittest.TestSuite()
+ loader = unittest.defaultTestLoader.loadTestsFromTestCase
+ testsuite.addTest(loader(TestKahan))
+ return testsuite
+
+
+if __name__ == '__main__':
+ runner = unittest.TextTestRunner()
+ runner.run(suite())
diff --git a/silx/opencl/test/test_stats.py b/silx/opencl/test/test_stats.py
new file mode 100644
index 0000000..b5127c8
--- /dev/null
+++ b/silx/opencl/test/test_stats.py
@@ -0,0 +1,114 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Project: Sift implementation in Python + OpenCL
+# https://github.com/silx-kit/silx
+#
+# Permission is hereby granted, free of charge, to any person
+# obtaining a copy of this software and associated documentation
+# files (the "Software"), to deal in the Software without
+# restriction, including without limitation the rights to use,
+# copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following
+# conditions:
+#
+# The above copyright notice and this permission notice shall be
+# included in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+# OTHER DEALINGS IN THE SOFTWARE.
+
+"""
+Simple test of an addition
+"""
+
+from __future__ import division, print_function
+
+__authors__ = ["Henri Payno, Jérôme Kieffer"]
+__contact__ = "jerome.kieffer@esrf.eu"
+__license__ = "MIT"
+__copyright__ = "2013 European Synchrotron Radiation Facility, Grenoble, France"
+__date__ = "13/12/2018"
+
+import logging
+import time
+import numpy
+
+import unittest
+from ..common import ocl
+if ocl:
+ import pyopencl
+ import pyopencl.array
+ from ..statistics import StatResults, Statistics
+from ..utils import get_opencl_code
+logger = logging.getLogger(__name__)
+
+
+@unittest.skipUnless(ocl, "PyOpenCl is missing")
+class TestStatistics(unittest.TestCase):
+
+ @classmethod
+ def setUpClass(cls):
+ cls.size = 1 << 20 # 1 million elements
+ cls.data = numpy.random.randint(0, 65000, cls.size).astype("uint16")
+ t0 = time.time()
+ cls.ref = StatResults(cls.data.min(), cls.data.max(), cls.data.size,
+ cls.data.sum(), cls.data.mean(), cls.data.std() ** 2,
+ cls.data.std())
+ t1 = time.time()
+ cls.ref_time = t1 - t0
+
+ @classmethod
+ def tearDownClass(cls):
+ cls.size = cls.ref = cls.data = cls.ref_time = None
+
+ @classmethod
+ def validate(cls, res):
+ return (
+ (res.min == cls.ref.min) and
+ (res.max == cls.ref.max) and
+ (res.cnt == cls.ref.cnt) and
+ abs(res.mean - cls.ref.mean) < 0.01 and
+ abs(res.std - cls.ref.std) < 0.1)
+
+ def test_measurement(self):
+ """
+ tests that all devices are working properly ...
+ """
+ logger.info("Reference results: %s", self.ref)
+ for pid, platform in enumerate(ocl.platforms):
+ for did, device in enumerate(platform.devices):
+ try:
+ s = Statistics(template=self.data, platformid=pid, deviceid=did)
+ except Exception as err:
+ failed_init = True
+ res = StatResults(0,0,0,0,0,0,0)
+ else:
+ failed_init = False
+ t0 = time.time()
+ res = s(self.data)
+ t1 = time.time()
+ logger.warning("failed_init %s", failed_init)
+ if failed_init or not self.validate(res):
+ logger.error("Failed on platform %s device %s", platform, device)
+ logger.error("Reference results: %s", self.ref)
+ logger.error("Faulty results: %s", res)
+ self.assertTrue(False, "Stat calculation failed on %s %s" % (platform, device))
+ logger.info("Runtime on %s/%s : %.3fms x%.1f", platform, device, 1000 * (t1 - t0), self.ref_time / (t1 - t0))
+
+
+def suite():
+ testSuite = unittest.TestSuite()
+ testSuite.addTest(TestStatistics("test_measurement"))
+ return testSuite
+
+
+if __name__ == '__main__':
+ unittest.main(defaultTest="suite")
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-04 19:33:20 +0000