diff options
Diffstat (limited to 'src/silx/math/test')
| -rw-r--r-- | src/silx/math/test/__init__.py | 23 | ||||
| -rw-r--r-- | src/silx/math/test/benchmark_combo.py | 192 | ||||
| -rw-r--r-- | src/silx/math/test/histo_benchmarks.py | 269 | ||||
| -rw-r--r-- | src/silx/math/test/test_HistogramndLut_nominal.py | 571 | ||||
| -rw-r--r-- | src/silx/math/test/test_calibration.py | 145 | ||||
| -rw-r--r-- | src/silx/math/test/test_colormap.py | 269 | ||||
| -rw-r--r-- | src/silx/math/test/test_combo.py | 207 | ||||
| -rw-r--r-- | src/silx/math/test/test_histogramnd_error.py | 519 | ||||
| -rw-r--r-- | src/silx/math/test/test_histogramnd_nominal.py | 937 | ||||
| -rw-r--r-- | src/silx/math/test/test_histogramnd_vs_np.py | 826 | ||||
| -rw-r--r-- | src/silx/math/test/test_interpolate.py | 125 | ||||
| -rw-r--r-- | src/silx/math/test/test_marchingcubes.py | 174 |
12 files changed, 4257 insertions, 0 deletions
diff --git a/src/silx/math/test/__init__.py b/src/silx/math/test/__init__.py new file mode 100644 index 0000000..ad9836c --- /dev/null +++ b/src/silx/math/test/__init__.py @@ -0,0 +1,23 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2019 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ diff --git a/src/silx/math/test/benchmark_combo.py b/src/silx/math/test/benchmark_combo.py new file mode 100644 index 0000000..c12f590 --- /dev/null +++ b/src/silx/math/test/benchmark_combo.py @@ -0,0 +1,192 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2017 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +"""Benchmarks of the combo module""" + +from __future__ import division + +__authors__ = ["T. Vincent"] +__license__ = "MIT" +__date__ = "17/01/2018" + + +import logging +import os.path +import time +import unittest + +import numpy + +from silx.test.utils import temp_dir +from silx.utils.testutils import ParametricTestCase + +from silx.math import combo + +_logger = logging.getLogger(__name__) +_logger.setLevel(logging.DEBUG) + + +class TestBenchmarkMinMax(ParametricTestCase): + """Benchmark of min max combo""" + + DTYPES = ('float32', 'float64', + 'int8', 'int16', 'int32', 'int64', + 'uint8', 'uint16', 'uint32', 'uint64') + + ARANGE = 'ascent', 'descent', 'random' + + EXPONENT = 3, 4, 5, 6, 7 + + def test_benchmark_min_max(self): + """Benchmark min_max without min positive. + + Compares with: + + - numpy.nanmin, numpy.nanmax and + - numpy.argmin, numpy.argmax + + It runs bench for different types, different data size and 3 + data sets: increasing , decreasing and random data. + """ + durations = {'min/max': [], 'argmin/max': [], 'combo': []} + + _logger.info('Benchmark against argmin/argmax and nanmin/nanmax') + + for dtype in self.DTYPES: + for arange in self.ARANGE: + for exponent in self.EXPONENT: + size = 10**exponent + with self.subTest(dtype=dtype, size=size, arange=arange): + if arange == 'ascent': + data = numpy.arange(0, size, 1, dtype=dtype) + elif arange == 'descent': + data = numpy.arange(size, 0, -1, dtype=dtype) + else: + if dtype in ('float32', 'float64'): + data = numpy.random.random(size) + else: + data = numpy.random.randint(10**6, size=size) + data = numpy.array(data, dtype=dtype) + + start = time.time() + ref_min = numpy.nanmin(data) + ref_max = numpy.nanmax(data) + durations['min/max'].append(time.time() - start) + + start = time.time() + ref_argmin = numpy.argmin(data) + ref_argmax = numpy.argmax(data) + durations['argmin/max'].append(time.time() - start) + + start = time.time() + result = combo.min_max(data, min_positive=False) + durations['combo'].append(time.time() - start) + + _logger.info( + '%s-%s-10**%d\tx%.2f argmin/max x%.2f min/max', + dtype, arange, exponent, + durations['argmin/max'][-1] / durations['combo'][-1], + durations['min/max'][-1] / durations['combo'][-1]) + + self.assertEqual(result.minimum, ref_min) + self.assertEqual(result.maximum, ref_max) + self.assertEqual(result.argmin, ref_argmin) + self.assertEqual(result.argmax, ref_argmax) + + self.show_results('min/max', durations, 'combo') + + def test_benchmark_min_pos(self): + """Benchmark min_max wit min positive. + + Compares with: + + - numpy.nanmin(data[data > 0]); numpy.nanmin(pos); numpy.nanmax(pos) + + It runs bench for different types, different data size and 3 + data sets: increasing , decreasing and random data. + """ + durations = {'min/max': [], 'combo': []} + + _logger.info('Benchmark against min, max, positive min') + + for dtype in self.DTYPES: + for arange in self.ARANGE: + for exponent in self.EXPONENT: + size = 10**exponent + with self.subTest(dtype=dtype, size=size, arange=arange): + if arange == 'ascent': + data = numpy.arange(0, size, 1, dtype=dtype) + elif arange == 'descent': + data = numpy.arange(size, 0, -1, dtype=dtype) + else: + if dtype in ('float32', 'float64'): + data = numpy.random.random(size) + else: + data = numpy.random.randint(10**6, size=size) + data = numpy.array(data, dtype=dtype) + + start = time.time() + ref_min_positive = numpy.nanmin(data[data > 0]) + ref_min = numpy.nanmin(data) + ref_max = numpy.nanmax(data) + durations['min/max'].append(time.time() - start) + + start = time.time() + result = combo.min_max(data, min_positive=True) + durations['combo'].append(time.time() - start) + + _logger.info( + '%s-%s-10**%d\tx%.2f min/minpos/max', + dtype, arange, exponent, + durations['min/max'][-1] / durations['combo'][-1]) + + self.assertEqual(result.min_positive, ref_min_positive) + self.assertEqual(result.minimum, ref_min) + self.assertEqual(result.maximum, ref_max) + + self.show_results('min/max/min positive', durations, 'combo') + + def show_results(self, title, durations, ref_key): + try: + from matplotlib import pyplot + except ImportError: + _logger.warning('matplotlib not available') + return + + pyplot.title(title) + pyplot.xlabel('-'.join(self.DTYPES)) + pyplot.ylabel('duration (sec)') + for label, values in durations.items(): + pyplot.semilogy(values, label=label) + pyplot.legend() + pyplot.show() + + pyplot.title(title) + pyplot.xlabel('-'.join(self.DTYPES)) + pyplot.ylabel('Duration ratio') + ref = numpy.array(durations[ref_key]) + for label, values in durations.items(): + values = numpy.array(values) + pyplot.plot(values/ref, label=label + ' / ' + ref_key) + pyplot.legend() + pyplot.show() diff --git a/src/silx/math/test/histo_benchmarks.py b/src/silx/math/test/histo_benchmarks.py new file mode 100644 index 0000000..7d3216d --- /dev/null +++ b/src/silx/math/test/histo_benchmarks.py @@ -0,0 +1,269 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +""" +histogramnd benchmarks, vs numpy.histogramdd (bin counts and weights). +""" + +import numpy as np + +import time + +from silx.math import histogramnd + + +def print_times(t0s, t1s, t2s, t3s): + c_times = t1s - t0s + np_times = t2s - t1s + np_w_times = t3s - t2s + + time_txt = 'min : {0: <7.3f}; max : {1: <7.3f}; avg : {2: <7.3f}' + + print('\tTimes :') + print('\tC : ' + time_txt.format(c_times.min(), + c_times.max(), + c_times.mean())) + print('\tNP : ' + time_txt.format(np_times.min(), + np_times.max(), + np_times.mean())) + print('\tNP(W) : ' + time_txt.format(np_w_times.min(), + np_w_times.max(), + np_w_times.mean())) + + +def commpare_results(txt, + times, + result_c, + result_np, + result_np_w, + sample, + weights, + raise_ex=False): + + if result_np: + hits_cmp = np.array_equal(result_c[0], result_np[0]) + else: + hits_cmp = None + + if result_np_w and result_c[1] is not None: + weights_cmp = np.array_equal(result_c[1], result_np_w[0]) + else: + weights_cmp = None + + if((hits_cmp is not None and not hits_cmp) or + (weights_cmp is not None and not weights_cmp)): + err_txt = (txt + ' : results arent the same : ' + 'hits : {0}, ' + 'weights : {1}.' + ''.format('OK' if hits_cmp else 'NOK', + 'OK' if weights_cmp else 'NOK')) + print('\t' + err_txt) + if raise_ex: + raise ValueError(err_txt) + return False + + result_txt = ' : results OK. c : {0: <7.3f};'.format(times[0]) + if result_np or result_np_w: + result_txt += (' np : {0: <7.3f}; ' + 'np (weights) {1: <7.3f}.' + ''.format(times[1], times[2])) + print('\t' + txt + result_txt) + return True + + +def benchmark(n_loops, + sample_shape, + sample_rng, + weights_rng, + histo_range, + n_bins, + weight_min, + weight_max, + last_bin_closed, + dtype=np.double, + do_weights=True, + do_numpy=True): + + int_min = 0 + int_max = 100000 + + sample = np.random.randint(int_min, + high=int_max, + size=sample_shape).astype(np.double) + sample = (sample_rng[0] + + (sample - int_min) * + (sample_rng[1] - sample_rng[0]) / + (int_max - int_min)) + sample = sample.astype(dtype) + + if do_weights: + weights = np.random.randint(int_min, + high=int_max, + size=(ssetup.pyample_shape[0],)) + weights = weights.astype(np.double) + weights = (weights_rng[0] + + (weights - int_min) * + (weights_rng[1] - weights_rng[0]) / + (int_max - int_min)) + else: + weights = None + + t0s = [] + t1s = [] + t2s = [] + t3s = [] + + for i in range(n_loops): + t0s.append(time.time()) + result_c = histogramnd(sample, + histo_range, + n_bins, + weights=weights, + weight_min=weight_min, + weight_max=weight_max, + last_bin_closed=last_bin_closed) + t1s.append(time.time()) + if do_numpy: + result_np = np.histogramdd(sample, + bins=n_bins, + range=histo_range) + t2s.append(time.time()) + result_np_w = np.histogramdd(sample, + bins=n_bins, + range=histo_range, + weights=weights) + t3s.append(time.time()) + else: + result_np = None + result_np_w = None + t2s.append(0) + t3s.append(0) + + commpare_results('Run {0}'.format(i), + [t1s[-1] - t0s[-1], t2s[-1] - t1s[-1], t3s[-1] - t2s[-1]], + result_c, + result_np, + result_np_w, + sample, + weights) + + print_times(np.array(t0s), np.array(t1s), np.array(t2s), np.array(t3s)) + + +def run_benchmark(dtype=np.double, + do_weights=True, + do_numpy=True): + n_loops = 5 + + weights_rng = [0., 100.] + sample_rng = [0., 100.] + + weight_min = None + weight_max = None + last_bin_closed = True + + # ==================================================== + # ==================================================== + # 1D + # ==================================================== + # ==================================================== + + print('==========================') + print(' 1D [{0}]'.format(dtype)) + print('==========================') + sample_shape = (10**7,) + histo_range = [[0., 100.]] + n_bins = 30 + + benchmark(n_loops, + sample_shape, + sample_rng, + weights_rng, + histo_range, + n_bins, + weight_min, + weight_max, + last_bin_closed, + dtype=dtype, + do_weights=True, + do_numpy=do_numpy) + + # ==================================================== + # ==================================================== + # 2D + # ==================================================== + # ==================================================== + + print('==========================') + print(' 2D [{0}]'.format(dtype)) + print('==========================') + sample_shape = (10**7, 2) + histo_range = [[0., 100.], [0., 100.]] + n_bins = 30 + + benchmark(n_loops, + sample_shape, + sample_rng, + weights_rng, + histo_range, + n_bins, + weight_min, + weight_max, + last_bin_closed, + dtype=dtype, + do_weights=True, + do_numpy=do_numpy) + + # ==================================================== + # ==================================================== + # 3D + # ==================================================== + # ==================================================== + + print('==========================') + print(' 3D [{0}]'.format(dtype)) + print('==========================') + sample_shape = (10**7, 3) + histo_range = np.array([[0., 100.], [0., 100.], [0., 100.]]) + n_bins = 30 + + benchmark(n_loops, + sample_shape, + sample_rng, + weights_rng, + histo_range, + n_bins, + weight_min, + weight_max, + last_bin_closed, + dtype=dtype, + do_weights=True, + do_numpy=do_numpy) + +if __name__ == '__main__': + types = (np.double, np.int32, np.float32,) + + for t in types: + run_benchmark(t, + do_weights=True, + do_numpy=True) diff --git a/src/silx/math/test/test_HistogramndLut_nominal.py b/src/silx/math/test/test_HistogramndLut_nominal.py new file mode 100644 index 0000000..52e003c --- /dev/null +++ b/src/silx/math/test/test_HistogramndLut_nominal.py @@ -0,0 +1,571 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2019 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +""" +Nominal tests of the HistogramndLut function. +""" + +import unittest + +import numpy as np + +from silx.math import HistogramndLut + + +def _get_bin_edges(histo_range, n_bins, n_dims): + edges = [] + for i_dim in range(n_dims): + edges.append(histo_range[i_dim, 0] + + np.arange(n_bins[i_dim] + 1) * + (histo_range[i_dim, 1] - histo_range[i_dim, 0]) / + n_bins[i_dim]) + return tuple(edges) + + +# ============================================================== +# ============================================================== +# ============================================================== + + +class _TestHistogramndLut_nominal(unittest.TestCase): + """ + Unit tests of the HistogramndLut class. + """ + __test__ = False # ignore abstract class + + ndims = None + + def setUp(self): + ndims = self.ndims + if ndims is None: + self.skipTest("Abstract class") + self.tested_dim = ndims-1 + + if ndims is None: + raise ValueError('ndims class member not set.') + + sample = np.array([5.5, -3.3, + 0., -0.5, + 3.3, 8.8, + -7.7, 6.0, + -4.0]) + + weights = np.array([500.5, -300.3, + 0.01, -0.5, + 300.3, 800.8, + -700.7, 600.6, + -400.4]) + + n_elems = len(sample) + + if ndims == 1: + shape = (n_elems,) + else: + shape = (n_elems, ndims) + + self.sample = np.zeros(shape=shape, dtype=sample.dtype) + if ndims == 1: + self.sample = sample + else: + self.sample[..., ndims-1] = sample + + self.weights = weights + + # the tests are performed along one dimension, + # all the other bins indices along the other dimensions + # are expected to be 2 + # (e.g : when testing a 2D sample : [0, x] will go into + # bin [2, y] because of the bin ranges [-2, 2] and n_bins = 4 + # for the first dimension) + self.other_axes_index = 2 + self.histo_range = np.repeat([[-2., 2.]], ndims, axis=0) + self.histo_range[ndims-1] = [-4., 6.] + + self.n_bins = np.array([4]*ndims) + self.n_bins[ndims-1] = 5 + + if ndims == 1: + def fill_histo(h, v, dim, op=None): + if op: + h[:] = op(h[:], v) + else: + h[:] = v + self.fill_histo = fill_histo + else: + def fill_histo(h, v, dim, op=None): + idx = [self.other_axes_index]*len(h.shape) + idx[dim] = slice(0, None) + idx = tuple(idx) + if op: + h[idx] = op(h[idx], v) + else: + h[idx] = v + self.fill_histo = fill_histo + + def test_nominal_bin_edges(self): + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + bin_edges = instance.bins_edges + + expected_edges = _get_bin_edges(self.histo_range, + self.n_bins, + self.ndims) + + for i_edges, edges in enumerate(expected_edges): + self.assertTrue(np.array_equal(bin_edges[i_edges], + expected_edges[i_edges]), + msg='Testing bin_edges for dim {0}' + ''.format(i_edges+1)) + + def test_nominal_histo_range(self): + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + histo_range = instance.histo_range + + self.assertTrue(np.array_equal(histo_range, self.histo_range)) + + def test_nominal_last_bin_closed(self): + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + last_bin_closed = instance.last_bin_closed + + self.assertEqual(last_bin_closed, False) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins, + last_bin_closed=True) + + last_bin_closed = instance.last_bin_closed + + self.assertEqual(last_bin_closed, True) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins, + last_bin_closed=False) + + last_bin_closed = instance.last_bin_closed + + self.assertEqual(last_bin_closed, False) + + def test_nominal_n_bins_array(self): + + test_n_bins = np.arange(self.ndims) + 10 + instance = HistogramndLut(self.sample, + self.histo_range, + test_n_bins) + + n_bins = instance.n_bins + + self.assertTrue(np.array_equal(test_n_bins, n_bins)) + + def test_nominal_n_bins_scalar(self): + + test_n_bins = 10 + expected_n_bins = np.array([test_n_bins] * self.ndims) + instance = HistogramndLut(self.sample, + self.histo_range, + test_n_bins) + + n_bins = instance.n_bins + + self.assertTrue(np.array_equal(expected_n_bins, n_bins)) + + def test_nominal_histo_ref(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + instance.accumulate(self.weights) + + histo = instance.histo() + w_histo = instance.weighted_histo() + histo_ref = instance.histo(copy=False) + w_histo_ref = instance.weighted_histo(copy=False) + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + self.assertTrue(np.array_equal(histo_ref, expected_h)) + self.assertTrue(np.array_equal(w_histo_ref, expected_c)) + + histo_ref[0, ...] = histo_ref[0, ...] + 10 + w_histo_ref[0, ...] = w_histo_ref[0, ...] + 20 + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + self.assertFalse(np.array_equal(histo_ref, expected_h)) + self.assertFalse(np.array_equal(w_histo_ref, expected_c)) + + histo_2 = instance.histo() + w_histo_2 = instance.weighted_histo() + + self.assertFalse(np.array_equal(histo_2, expected_h)) + self.assertFalse(np.array_equal(w_histo_2, expected_c)) + self.assertTrue(np.array_equal(histo_2, histo_ref)) + self.assertTrue(np.array_equal(w_histo_2, w_histo_ref)) + + def test_nominal_accumulate_once(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + instance.accumulate(self.weights) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + self.assertEqual(w_histo.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + self.assertTrue(np.array_equal(instance.histo(), expected_h)) + self.assertTrue(np.array_equal(instance.weighted_histo(), + expected_c)) + + def test_nominal_accumulate_twice(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + # calling accumulate twice + expected_h *= 2 + expected_c *= 2 + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + instance.accumulate(self.weights) + + instance.accumulate(self.weights) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + self.assertEqual(w_histo.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + self.assertTrue(np.array_equal(instance.histo(), expected_h)) + self.assertTrue(np.array_equal(instance.weighted_histo(), + expected_c)) + + def test_nominal_apply_lut_once(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + histo, w_histo = instance.apply_lut(self.weights) + + self.assertEqual(w_histo.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + self.assertEqual(instance.histo(), None) + self.assertEqual(instance.weighted_histo(), None) + + def test_nominal_apply_lut_twice(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + # calling apply_lut twice + expected_h *= 2 + expected_c *= 2 + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + histo, w_histo = instance.apply_lut(self.weights) + histo_2, w_histo_2 = instance.apply_lut(self.weights, + histo=histo, + weighted_histo=w_histo) + + self.assertEqual(id(histo), id(histo_2)) + self.assertEqual(id(w_histo), id(w_histo_2)) + self.assertEqual(w_histo.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + self.assertEqual(instance.histo(), None) + self.assertEqual(instance.weighted_histo(), None) + + def test_nominal_accumulate_last_bin_closed(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 2]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 1101.1]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins, + last_bin_closed=True) + + instance.accumulate(self.weights) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + self.assertEqual(w_histo.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + + def test_nominal_accumulate_weight_min_max(self): + """ + """ + weight_min = -299.9 + weight_max = 499.9 + + expected_h_tpl = np.array([0, 1, 1, 1, 0]) + expected_c_tpl = np.array([0., -0.5, 0.01, 300.3, 0.]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + instance.accumulate(self.weights, + weight_min=weight_min, + weight_max=weight_max) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + self.assertEqual(w_histo.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + + def test_nominal_accumulate_forced_int32(self): + """ + double weights, int32 weighted_histogram + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700, 0, 0, 300, 500]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins, + dtype=np.int32) + + instance.accumulate(self.weights) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + self.assertEqual(w_histo.dtype, np.int32) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + + def test_nominal_accumulate_forced_float32(self): + """ + int32 weights, float32 weighted_histogram + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700., 0., 0., 300., 500.]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.float32) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins, + dtype=np.float32) + + instance.accumulate(self.weights.astype(np.int32)) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + self.assertEqual(w_histo.dtype, np.float32) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + + def test_nominal_accumulate_int32(self): + """ + int32 weights + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700, 0, 0, 300, 500]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.int32) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + instance.accumulate(self.weights.astype(np.int32)) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + self.assertEqual(w_histo.dtype, np.int32) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + + def test_nominal_accumulate_int32_double(self): + """ + int32 weights + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700, 0, 0, 300, 500]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.int32) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + instance = HistogramndLut(self.sample, + self.histo_range, + self.n_bins) + + instance.accumulate(self.weights.astype(np.int32)) + instance.accumulate(self.weights) + + histo = instance.histo() + w_histo = instance.weighted_histo() + + expected_h *= 2 + expected_c *= 2 + + self.assertEqual(w_histo.dtype, np.int32) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(w_histo, expected_c)) + + def testNoneNativeTypes(self): + type = self.sample.dtype.newbyteorder("B") + sampleB = self.sample.astype(type) + + type = self.sample.dtype.newbyteorder("L") + sampleL = self.sample.astype(type) + + histo_inst = HistogramndLut(sampleB, + self.histo_range, + self.n_bins) + + histo_inst = HistogramndLut(sampleL, + self.histo_range, + self.n_bins) + + +class TestHistogramndLut_nominal_1d(_TestHistogramndLut_nominal): + __test__ = True # because _TestHistogramndLut_nominal is ignored + ndims = 1 + + +class TestHistogramndLut_nominal_2d(_TestHistogramndLut_nominal): + __test__ = True # because _TestHistogramndLut_nominal is ignored + ndims = 2 + + +class TestHistogramndLut_nominal_3d(_TestHistogramndLut_nominal): + __test__ = True # because _TestHistogramndLut_nominal is ignored + ndims = 3 diff --git a/src/silx/math/test/test_calibration.py b/src/silx/math/test/test_calibration.py new file mode 100644 index 0000000..7158293 --- /dev/null +++ b/src/silx/math/test/test_calibration.py @@ -0,0 +1,145 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2018 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +"""Tests of the calibration module""" + +from __future__ import division + +__authors__ = ["P. Knobel"] +__license__ = "MIT" +__date__ = "14/05/2018" + + +import unittest + +import numpy + +from silx.math.calibration import NoCalibration, LinearCalibration, \ + ArrayCalibration, FunctionCalibration + + +X = numpy.array([3.14, 2.73, 1337]) + + +class TestNoCalibration(unittest.TestCase): + def setUp(self): + self.calib = NoCalibration() + + def testIsAffine(self): + self.assertTrue(self.calib.is_affine()) + + def testSlope(self): + self.assertEqual(self.calib.get_slope(), 1.) + + def testYIntercept(self): + self.assertEqual(self.calib(0.), + 0.) + + def testCall(self): + self.assertTrue(numpy.array_equal(self.calib(X), X)) + + +class TestLinearCalibration(unittest.TestCase): + def setUp(self): + self.y_intercept = 1.5 + self.slope = 2.5 + self.calib = LinearCalibration(y_intercept=self.y_intercept, + slope=self.slope) + + def testIsAffine(self): + self.assertTrue(self.calib.is_affine()) + + def testSlope(self): + self.assertEqual(self.calib.get_slope(), self.slope) + + def testYIntercept(self): + self.assertEqual(self.calib(0.), + self.y_intercept) + + def testCall(self): + self.assertTrue(numpy.array_equal(self.calib(X), + self.y_intercept + self.slope * X)) + + +class TestArrayCalibration(unittest.TestCase): + def setUp(self): + self.arr = numpy.array([45.2, 25.3, 666., -8.]) + self.calib = ArrayCalibration(self.arr) + self.affine_calib = ArrayCalibration([0.1, 0.2, 0.3]) + + def testIsAffine(self): + self.assertFalse(self.calib.is_affine()) + self.assertTrue(self.affine_calib.is_affine()) + + def testSlope(self): + with self.assertRaises(AttributeError): + self.calib.get_slope() + self.assertEqual(self.affine_calib.get_slope(), + 0.1) + + def testYIntercept(self): + self.assertEqual(self.calib(0), + self.arr[0]) + + def testCall(self): + with self.assertRaises(ValueError): + # X is an array with a different shape + self.calib(X) + + with self.assertRaises(ValueError): + # floats are not valid indices + self.calib(3.14) + + self.assertTrue( + numpy.array_equal(self.calib([1, 2, 3, 4]), + self.arr)) + + for idx, value in enumerate(self.arr): + self.assertEqual(self.calib(idx), value) + + +class TestFunctionCalibration(unittest.TestCase): + def setUp(self): + self.non_affine_fun = numpy.sin + self.non_affine_calib = FunctionCalibration(self.non_affine_fun) + + self.affine_fun = lambda x: 52. * x + 0.01 + self.affine_calib = FunctionCalibration(self.affine_fun, + is_affine=True) + + def testIsAffine(self): + self.assertFalse(self.non_affine_calib.is_affine()) + self.assertTrue(self.affine_calib.is_affine()) + + def testSlope(self): + with self.assertRaises(AttributeError): + self.non_affine_calib.get_slope() + self.assertAlmostEqual(self.affine_calib.get_slope(), + 52.) + + def testCall(self): + for x in X: + self.assertAlmostEqual(self.non_affine_calib(x), + self.non_affine_fun(x)) + self.assertAlmostEqual(self.affine_calib(x), + self.affine_fun(x)) diff --git a/src/silx/math/test/test_colormap.py b/src/silx/math/test/test_colormap.py new file mode 100644 index 0000000..0b0ec59 --- /dev/null +++ b/src/silx/math/test/test_colormap.py @@ -0,0 +1,269 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2018-2021 European Synchrotron Radiation Facility +# +# 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 for colormap mapping implementation""" + +from __future__ import division + +__authors__ = ["T. Vincent"] +__license__ = "MIT" +__date__ = "16/05/2018" + + +import logging +import sys + +import numpy + +from silx.utils.testutils import ParametricTestCase +from silx.math import colormap + + +_logger = logging.getLogger(__name__) + + +class TestNormalization(ParametricTestCase): + """Test silx.math.colormap.Normalization sub classes""" + + def _testCodec(self, normalization, rtol=1e-5): + """Test apply/revert for normalizations""" + test_data = (numpy.arange(1, 10, dtype=numpy.int32), + numpy.linspace(1., 100., 1000, dtype=numpy.float32), + numpy.linspace(-1., 1., 100, dtype=numpy.float32), + 1., + 1) + + for index in range(len(test_data)): + with self.subTest(normalization=normalization, data_index=index): + data = test_data[index] + normalized = normalization.apply(data, 1., 100.) + result = normalization.revert(normalized, 1., 100.) + + self.assertTrue(numpy.array_equal( + numpy.isnan(normalized), numpy.isnan(result))) + + if isinstance(data, numpy.ndarray): + notNaN = numpy.logical_not(numpy.isnan(result)) + data = data[notNaN] + result = result[notNaN] + self.assertTrue(numpy.allclose(data, result, rtol=rtol)) + + def testLinearNormalization(self): + """Test for LinearNormalization""" + normalization = colormap.LinearNormalization() + self._testCodec(normalization) + + def testLogarithmicNormalization(self): + """Test for LogarithmicNormalization""" + normalization = colormap.LogarithmicNormalization() + # relative tolerance is higher because of the log approximation + self._testCodec(normalization, rtol=1e-3) + + # Specific extra tests + self.assertTrue(numpy.isnan(normalization.apply(-1., 1., 100.))) + self.assertTrue(numpy.isnan(normalization.apply(numpy.nan, 1., 100.))) + self.assertEqual(normalization.apply(numpy.inf, 1., 100.), numpy.inf) + self.assertEqual(normalization.apply(0, 1., 100.), - numpy.inf) + + def testArcsinhNormalization(self): + """Test for ArcsinhNormalization""" + self._testCodec(colormap.ArcsinhNormalization()) + + def testSqrtNormalization(self): + """Test for SqrtNormalization""" + normalization = colormap.SqrtNormalization() + self._testCodec(normalization) + + # Specific extra tests + self.assertTrue(numpy.isnan(normalization.apply(-1., 0., 100.))) + self.assertTrue(numpy.isnan(normalization.apply(numpy.nan, 0., 100.))) + self.assertEqual(normalization.apply(numpy.inf, 0., 100.), numpy.inf) + self.assertEqual(normalization.apply(0, 0., 100.), 0.) + + +class TestColormap(ParametricTestCase): + """Test silx.math.colormap.cmap""" + + NORMALIZATIONS = ( + 'linear', + 'log', + 'arcsinh', + 'sqrt', + colormap.LinearNormalization(), + colormap.LogarithmicNormalization(), + colormap.GammaNormalization(2.), + colormap.GammaNormalization(0.5)) + + @staticmethod + def ref_colormap(data, colors, vmin, vmax, normalization, nan_color): + """Reference implementation of colormap + + :param numpy.ndarray data: Data to convert + :param numpy.ndarray colors: Color look-up-table + :param float vmin: Lower bound of the colormap range + :param float vmax: Upper bound of the colormap range + :param str normalization: Normalization to use + :param Union[numpy.ndarray, None] nan_color: Color to use for NaN + """ + norm_functions = {'linear': lambda v: v, + 'log': numpy.log10, + 'arcsinh': numpy.arcsinh, + 'sqrt': numpy.sqrt} + + if isinstance(normalization, str): + norm_function = norm_functions[normalization] + else: + def norm_function(value): + return normalization.apply(value, vmin, vmax) + + with numpy.errstate(divide='ignore', invalid='ignore'): + # Ignore divide by zero and invalid value encountered in log10, sqrt + norm_data, vmin, vmax = map(norm_function, (data, vmin, vmax)) + + if normalization == 'arcsinh' and sys.platform == 'win32': + # There is a difference of behavior of numpy.arcsinh + # between Windows and other OS for results of infinite values + # This makes Windows behaves as Linux and MacOS + norm_data[data == numpy.inf] = numpy.inf + norm_data[data == -numpy.inf] = -numpy.inf + + nb_colors = len(colors) + scale = nb_colors / (vmax - vmin) + + # Substraction must be done in float to avoid overflow with uint + indices = numpy.clip(scale * (norm_data - float(vmin)), + 0, nb_colors - 1) + indices[numpy.isnan(indices)] = nb_colors # Use an extra index for NaN + indices = indices.astype('uint') + + # Add NaN color to array + if nan_color is None: + nan_color = (0,) * colors.shape[-1] + colors = numpy.append(colors, numpy.atleast_2d(nan_color), axis=0) + + return colors[indices] + + def _test(self, data, colors, vmin, vmax, normalization, nan_color): + """Run test of colormap against alternative implementation + + :param numpy.ndarray data: Data to convert + :param numpy.ndarray colors: Color look-up-table + :param float vmin: Lower bound of the colormap range + :param float vmax: Upper bound of the colormap range + :param str normalization: Normalization to use + :param Union[numpy.ndarray, None] nan_color: Color to use for NaN + """ + image = colormap.cmap( + data, colors, vmin, vmax, normalization, nan_color) + + ref_image = self.ref_colormap( + data, colors, vmin, vmax, normalization, nan_color) + + self.assertTrue(numpy.allclose(ref_image, image)) + self.assertEqual(image.dtype, colors.dtype) + self.assertEqual(image.shape, data.shape + (colors.shape[-1],)) + + def test(self): + """Test all dtypes with finite data + + Test all supported types and endianness + """ + colors = numpy.zeros((256, 4), dtype=numpy.uint8) + colors[:, 0] = numpy.arange(len(colors)) + colors[:, 3] = 255 + + # Generates (u)int and floats types + dtypes = [e + k + i for e in '<>' for k in 'uif' for i in '1248' + if k != 'f' or i != '1'] + dtypes.append(numpy.dtype(numpy.longdouble).name) # Add long double + + for normalization in self.NORMALIZATIONS: + for dtype in dtypes: + with self.subTest(dtype=dtype, normalization=normalization): + _logger.info('normalization: %s, dtype: %s', + normalization, dtype) + data = numpy.arange(-5, 15, dtype=dtype).reshape(4, 5) + + self._test(data, colors, 1, 10, normalization, None) + + def test_not_finite(self): + """Test float data with not finite values""" + colors = numpy.zeros((256, 4), dtype=numpy.uint8) + colors[:, 0] = numpy.arange(len(colors)) + colors[:, 3] = 255 + + test_data = { # message: data + 'no finite values': (float('inf'), float('-inf'), float('nan')), + 'only NaN': (float('nan'), float('nan'), float('nan')), + 'mix finite/not finite': (float('inf'), float('-inf'), 1., float('nan')), + } + + for normalization in self.NORMALIZATIONS: + for msg, data in test_data.items(): + with self.subTest(msg, normalization=normalization): + _logger.info('normalization: %s, %s', normalization, msg) + data = numpy.array(data, dtype=numpy.float64) + self._test(data, colors, 1, 10, normalization, (0, 0, 0, 0)) + + def test_errors(self): + """Test raising exception for bad vmin, vmax, normalization parameters + """ + colors = numpy.zeros((256, 4), dtype=numpy.uint8) + colors[:, 0] = numpy.arange(len(colors)) + colors[:, 3] = 255 + + data = numpy.arange(10, dtype=numpy.float64) + + test_params = [ # (vmin, vmax, normalization) + (-1., 2., 'log'), + (0., 1., 'log'), + (1., 0., 'log'), + (-1., 1., 'sqrt'), + (1., -1., 'sqrt'), + ] + + for vmin, vmax, normalization in test_params: + with self.subTest( + vmin=vmin, vmax=vmax, normalization=normalization): + _logger.info('normalization: %s, range: [%f, %f]', + normalization, vmin, vmax) + with self.assertRaises(ValueError): + self._test(data, colors, vmin, vmax, normalization, None) + + +def test_apply_colormap(): + """Basic test of silx.math.colormap.apply_colormap""" + data = numpy.arange(256) + expected_colors = numpy.empty((256, 4), dtype=numpy.uint8) + expected_colors[:, :3] = numpy.arange(256, dtype=numpy.uint8).reshape(256, 1) + expected_colors[:, 3] = 255 + colors = colormap.apply_colormap( + data, + colormap="gray", + norm="linear", + autoscale="minmax", + vmin=None, + vmax=None, + gamma=1.0) + assert numpy.array_equal(colors, expected_colors) diff --git a/src/silx/math/test/test_combo.py b/src/silx/math/test/test_combo.py new file mode 100644 index 0000000..9a96923 --- /dev/null +++ b/src/silx/math/test/test_combo.py @@ -0,0 +1,207 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2020 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +"""Tests of the combo module""" + +from __future__ import division + +__authors__ = ["T. Vincent"] +__license__ = "MIT" +__date__ = "17/01/2018" + + +import unittest + +import numpy + +from silx.utils.testutils import ParametricTestCase + +from silx.math.combo import min_max + + +class TestMinMax(ParametricTestCase): + """Tests of min max combo""" + + FLOATING_DTYPES = 'float32', 'float64' + if hasattr(numpy, "float128"): + FLOATING_DTYPES += ('float128',) + SIGNED_INT_DTYPES = 'int8', 'int16', 'int32', 'int64' + UNSIGNED_INT_DTYPES = 'uint8', 'uint16', 'uint32', 'uint64' + DTYPES = FLOATING_DTYPES + SIGNED_INT_DTYPES + UNSIGNED_INT_DTYPES + + def _numpy_min_max(self, data, min_positive=False, finite=False): + """Reference numpy implementation of min_max + + :param numpy.ndarray data: Data set to use for test + :param bool min_positive: True to test with positive min + :param bool finite: True to only test finite values + """ + data = numpy.array(data, copy=False) + if data.size == 0: + raise ValueError('Zero-sized array') + + minimum = None + argmin = None + maximum = None + argmax = None + min_pos = None + argmin_pos = None + + if finite: + filtered_data = data[numpy.isfinite(data)] + else: + filtered_data = data + + if filtered_data.size > 0: + if numpy.all(numpy.isnan(filtered_data)): + minimum = numpy.nan + argmin = 0 + maximum = numpy.nan + argmax = 0 + else: + minimum = numpy.nanmin(filtered_data) + # nanargmin equivalent + argmin = numpy.where(data == minimum)[0][0] + maximum = numpy.nanmax(filtered_data) + # nanargmax equivalent + argmax = numpy.where(data == maximum)[0][0] + + if min_positive: + with numpy.errstate(invalid='ignore'): + # Ignore invalid value encountered in greater + pos_data = filtered_data[filtered_data > 0] + if pos_data.size > 0: + min_pos = numpy.min(pos_data) + argmin_pos = numpy.where(data == min_pos)[0][0] + + return minimum, min_pos, maximum, argmin, argmin_pos, argmax + + def _test_min_max(self, data, min_positive, finite=False): + """Compare min_max with numpy for the given dataset + + :param numpy.ndarray data: Data set to use for test + :param bool min_positive: True to test with positive min + :param bool finite: True to only test finite values + """ + minimum, min_pos, maximum, argmin, argmin_pos, argmax = \ + self._numpy_min_max(data, min_positive, finite) + + result = min_max(data, min_positive, finite) + + self.assertSimilar(minimum, result.minimum) + self.assertSimilar(min_pos, result.min_positive) + self.assertSimilar(maximum, result.maximum) + self.assertSimilar(argmin, result.argmin) + self.assertSimilar(argmin_pos, result.argmin_positive) + self.assertSimilar(argmax, result.argmax) + + def assertSimilar(self, a, b): + """Assert that a and b are both None or NaN or that a == b.""" + self.assertTrue((a is None and b is None) or + (numpy.isnan(a) and numpy.isnan(b)) or + a == b) + + def test_different_datasets(self): + """Test min_max with different numpy.arange datasets.""" + size = 1000 + + for dtype in self.DTYPES: + + tests = { + '0 to N': (0, 1), + 'N-1 to 0': (size - 1, -1)} + if dtype not in self.UNSIGNED_INT_DTYPES: + tests['N/2 to -N/2'] = size // 2, -1 + tests['0 to -N'] = 0, -1 + + for name, (start, step) in tests.items(): + for min_positive in (True, False): + with self.subTest(dtype=dtype, + min_positive=min_positive, + data=name): + data = numpy.arange( + start, start + step * size, step, dtype=dtype) + + self._test_min_max(data, min_positive) + + def test_nodata(self): + """Test min_max with None and empty array""" + for dtype in self.DTYPES: + with self.subTest(dtype=dtype): + with self.assertRaises(TypeError): + min_max(None) + + data = numpy.array((), dtype=dtype) + with self.assertRaises(ValueError): + min_max(data) + + NAN_TEST_DATA = [ + (float('nan'), float('nan')), # All NaNs + (float('nan'), 1.0), # NaN first and positive + (float('nan'), -1.0), # NaN first and negative + (1.0, 2.0, float('nan')), # NaN last and positive + (-1.0, -2.0, float('nan')), # NaN last and negative + (1.0, float('nan'), -1.0), # Some NaN + ] + + def test_nandata(self): + """Test min_max with NaN in data""" + for dtype in self.FLOATING_DTYPES: + for data in self.NAN_TEST_DATA: + with self.subTest(dtype=dtype, data=data): + data = numpy.array(data, dtype=dtype) + self._test_min_max(data, min_positive=True) + + INF_TEST_DATA = [ + [float('inf')] * 3, # All +inf + [float('-inf')] * 3, # All -inf + (float('inf'), float('-inf')), # + and - inf + (float('inf'), float('-inf'), float('nan')), # +/-inf, nan last + (float('nan'), float('-inf'), float('inf')), # +/-inf, nan first + (float('inf'), float('nan'), float('-inf')), # +/-inf, nan center + ] + + def test_infdata(self): + """Test min_max with inf.""" + for dtype in self.FLOATING_DTYPES: + for data in self.INF_TEST_DATA: + with self.subTest(dtype=dtype, data=data): + data = numpy.array(data, dtype=dtype) + self._test_min_max(data, min_positive=True) + + def test_finite(self): + """Test min_max with finite=True""" + tests = [ + (-1., 2., 0.), # Basic test + (float('nan'), float('inf'), float('-inf')), # NaN + Inf + (float('nan'), float('inf'), -2, float('-inf')), # NaN + Inf + 1 value + (float('inf'), -3, -2), # values + inf + ] + tests += self.INF_TEST_DATA + tests += self.NAN_TEST_DATA + + for dtype in self.FLOATING_DTYPES: + for data in tests: + with self.subTest(dtype=dtype, data=data): + data = numpy.array(data, dtype=dtype) + self._test_min_max(data, min_positive=True, finite=True) diff --git a/src/silx/math/test/test_histogramnd_error.py b/src/silx/math/test/test_histogramnd_error.py new file mode 100644 index 0000000..22304cb --- /dev/null +++ b/src/silx/math/test/test_histogramnd_error.py @@ -0,0 +1,519 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ + +__authors__ = ["D. Naudet"] +__license__ = "MIT" +__date__ = "01/02/2016" + +""" +Tests of the histogramnd function, error cases. +""" +import sys +import platform +import unittest + +import numpy as np + +from silx.math.chistogramnd import chistogramnd as histogramnd +from silx.math import Histogramnd + + +# ============================================================== +# ============================================================== +# ============================================================== + + +class _Test_chistogramnd_errors(unittest.TestCase): + """ + Unit tests of the chistogramnd error cases. + """ + __test__ = False # ignore abstract class + + def setUp(self): + self.skipTest("Abstract class") + + def test_weights_shape(self): + """ + """ + + for err_w_shape in self.err_weights_shapes: + test_msg = ('Testing invalid weights shape : {0}' + ''.format(err_w_shape)) + + err_weights = np.random.randint(0, + high=10, + size=err_w_shape) + err_weights = err_weights.astype(np.double) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=err_weights)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, + '<weights> must be an array whose length ' + 'is equal to the number of samples.') + + def test_histo_range_shape(self): + """ + """ + n_dims = 1 if len(self.s_shape) == 1 else self.s_shape[1] + expected_txt_tpl = ('<histo_range> error : expected {n_dims} sets ' + 'of lower and upper bin edges, ' + 'got the following instead : {histo_range}. ' + '(provided <sample> contains ' + '{n_dims}D values)') + + for err_histo_range in self.err_histo_range_shapes: + test_msg = ('Testing invalid histo_range shape : {0}' + ''.format(err_histo_range)) + + expected_txt = expected_txt_tpl.format(histo_range=err_histo_range, + n_dims=n_dims) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + err_histo_range, + self.n_bins, + weights=self.weights)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_nbins_shape(self): + """ + """ + + expected_txt = ('n_bins must be either a scalar (same number ' + 'of bins for all dimensions) or ' + 'an array (number of bins for each ' + 'dimension).') + + for err_n_bins in self.err_n_bins_shapes: + test_msg = ('Testing invalid n_bins shape : {0}' + ''.format(err_n_bins)) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + err_n_bins, + weights=self.weights)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_nbins_values(self): + """ + """ + expected_txt = ('<n_bins> : only positive values allowed.') + + for err_n_bins in self.err_n_bins_values: + test_msg = ('Testing invalid n_bins value : {0}' + ''.format(err_n_bins)) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + err_n_bins, + weights=self.weights)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_histo_shape(self): + """ + """ + for err_h_shape in self.err_histo_shapes: + + # windows & python 2.7 : numpy shapes are long values + if platform.system() == 'Windows': + version = (sys.version_info.major, sys.version_info.minor) + if version <= (2, 7): + err_h_shape = tuple([long(val) for val in err_h_shape]) + + test_msg = ('Testing invalid histo shape : {0}' + ''.format(err_h_shape)) + + expected_txt = ('Provided <histo> array doesn\'t have ' + 'a shape compatible with <n_bins> ' + ': should be {0} instead of {1}.' + ''.format(self.h_shape, err_h_shape)) + + histo = np.zeros(shape=err_h_shape, dtype=np.uint32) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + histo=histo)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_histo_dtype(self): + """ + """ + for err_h_dtype in self.err_histo_dtypes: + test_msg = ('Testing invalid histo dtype : {0}' + ''.format(err_h_dtype)) + + histo = np.zeros(shape=self.h_shape, dtype=err_h_dtype) + + expected_txt = ('Provided <histo> array doesn\'t have ' + 'the expected type ' + ': should be {0} instead of {1}.' + ''.format(np.uint32, histo.dtype)) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + histo=histo)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_weighted_histo_shape(self): + """ + """ + # using the same values as histo + for err_h_shape in self.err_histo_shapes: + + # windows & python 2.7 : numpy shapes are long values + if platform.system() == 'Windows': + version = (sys.version_info.major, sys.version_info.minor) + if version <= (2, 7): + err_h_shape = tuple([long(val) for val in err_h_shape]) + + test_msg = ('Testing invalid weighted_histo shape : {0}' + ''.format(err_h_shape)) + + expected_txt = ('Provided <weighted_histo> array doesn\'t have ' + 'a shape compatible with <n_bins> ' + ': should be {0} instead of {1}.' + ''.format(self.h_shape, err_h_shape)) + + cumul = np.zeros(shape=err_h_shape, dtype=np.double) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + weighted_histo=cumul)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_cumul_dtype(self): + """ + """ + # using the same values as histo + for err_h_dtype in self.err_histo_dtypes: + test_msg = ('Testing invalid weighted_histo dtype : {0}' + ''.format(err_h_dtype)) + + cumul = np.zeros(shape=self.h_shape, dtype=err_h_dtype) + + expected_txt = ('Provided <weighted_histo> array doesn\'t have ' + 'the expected type ' + ': should be {0} or {1} instead of {2}.' + ''.format(np.float64, np.float32, cumul.dtype)) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + weighted_histo=cumul)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_wh_histo_dtype(self): + """ + """ + # using the same values as histo + for err_h_dtype in self.err_histo_dtypes: + test_msg = ('Testing invalid wh_dtype dtype : {0}' + ''.format(err_h_dtype)) + + expected_txt = ('<wh_dtype> type not supported : {0}.' + ''.format(err_h_dtype)) + + ex_str = None + try: + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + wh_dtype=err_h_dtype)[0:2] + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_unmanaged_dtypes(self): + """ + """ + for err_unmanaged_dtype in self.err_unmanaged_dtypes: + test_msg = ('Testing unmanaged dtypes : {0}' + ''.format(err_unmanaged_dtype)) + + sample = self.sample.astype(err_unmanaged_dtype[0]) + weights = self.weights.astype(err_unmanaged_dtype[1]) + + expected_txt = ('Case not supported - sample:{0} ' + 'and weights:{1}.' + ''.format(sample.dtype, + weights.dtype)) + + ex_str = None + try: + histogramnd(sample, + self.histo_range, + self.n_bins, + weights=weights) + except TypeError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str, msg=test_msg) + self.assertEqual(ex_str, expected_txt, msg=test_msg) + + def test_uncontiguous_histo(self): + """ + """ + # non contiguous array + shape = np.array(self.n_bins, ndmin=1) + shape[0] *= 2 + histo_tmp = np.zeros(shape) + histo = histo_tmp[::2, ...] + + expected_txt = ('<histo> must be a C_CONTIGUOUS numpy array.') + + ex_str = None + try: + histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + histo=histo) + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str) + self.assertEqual(ex_str, expected_txt) + + def test_uncontiguous_weighted_histo(self): + """ + """ + # non contiguous array + shape = np.array(self.n_bins, ndmin=1) + shape[0] *= 2 + cumul_tmp = np.zeros(shape) + cumul = cumul_tmp[::2, ...] + + expected_txt = ('<weighted_histo> must be a C_CONTIGUOUS numpy array.') + + ex_str = None + try: + histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + weighted_histo=cumul) + except ValueError as ex: + ex_str = str(ex) + + self.assertIsNotNone(ex_str) + self.assertEqual(ex_str, expected_txt) + + +class Test_chistogramnd_1D_errors(_Test_chistogramnd_errors): + """ + Unit tests of the 1D histogramnd error cases. + """ + __test__ = True # because _Test_chistogramnd_errors is ignored + + def setUp(self): + # nominal values + self.n_elements = 1000 + self.s_shape = (self.n_elements,) + self.w_shape = (self.n_elements,) + + self.histo_range = [0., 100.] + self.n_bins = 10 + + self.h_shape = (self.n_bins,) + + self.sample = np.random.randint(0, + high=10, + size=self.s_shape) + self.sample = self.sample.astype(np.double) + + self.weights = np.random.randint(0, + high=10, + size=self.w_shape) + self.weights = self.weights.astype(np.double) + + self.err_weights_shapes = ((self.n_elements+1,), + (self.n_elements-1,), + (self.n_elements-1, 3)) + self.err_histo_range_shapes = ([0.], + [0., 1., 2.], + [[0.], [1.]]) + self.err_n_bins_shapes = ([10, 2], + [[10], [2]]) + self.err_n_bins_values = (0, + [-10], + None) + self.err_histo_shapes = ((self.n_bins+1,), + (self.n_bins-1,), + (self.n_bins, self.n_bins)) + # these are used for testing the histo parameter as well + # as the weighted_histo parameter. + self.err_histo_dtypes = (np.uint16, + np.float16) + + self.err_unmanaged_dtypes = ((np.double, np.uint16), + (np.uint16, np.double), + (np.uint16, np.uint16)) + +class Test_chistogramnd_ND_range(unittest.TestCase): + """ + + """ + + def test_invalid_histo_range(self): + data = np.random.random((60, 60)) + nbins = 10 + + with self.assertRaises(ValueError): + histo_range = data.min(), np.inf + + Histogramnd(sample=data.ravel(), + histo_range=histo_range, + n_bins=nbins) + + histo_range = data.min(), np.nan + + Histogramnd(sample=data.ravel(), + histo_range=histo_range, + n_bins=nbins) + + +class Test_chistogramnd_ND_errors(_Test_chistogramnd_errors): + """ + Unit tests of the 3D histogramnd error cases. + """ + __test__ = True # because _Test_chistogramnd_errors is ignored + + def setUp(self): + # nominal values + self.n_elements = 1000 + self.s_shape = (self.n_elements, 3) + self.w_shape = (self.n_elements,) + + self.histo_range = [[0., 100.], [0., 100.], [0., 100.]] + self.n_bins = (10, 20, 30) + + self.h_shape = self.n_bins + + self.sample = np.random.randint(0, + high=10, + size=self.s_shape) + self.sample = self.sample.astype(np.double) + + self.weights = np.random.randint(0, + high=10, + size=self.w_shape) + self.weights = self.weights.astype(np.double) + + self.err_weights_shapes = ((self.n_elements+1,), + (self.n_elements-1,), + (self.n_elements-1, 3)) + self.err_histo_range_shapes = ([0.], + [0., 1.], + [[0., 10.], [0., 10.]], + [0., 10., 0, 10., 0, 10.]) + self.err_n_bins_shapes = ([10, 2], + [[10], [20], [30]]) + self.err_n_bins_values = (0, + [-10], + [10, 20, -4], + None, + [10, None, 30]) + self.err_histo_shapes = ((self.n_bins[0]+1, + self.n_bins[1], + self.n_bins[2]), + (self.n_bins[0], + self.n_bins[1], + self.n_bins[2]-1), + (self.n_bins[0], + self.n_bins[1]), + (self.n_bins[1], + self.n_bins[0], + self.n_bins[2]), + (self.n_bins[0], + self.n_bins[1], + self.n_bins[2], + 10) + ) + # these are used for testing the histo parameter as well + # as the weighted_histo parameter. + self.err_histo_dtypes = (np.uint16, + np.float16) + + self.err_unmanaged_dtypes = ((np.double, np.uint16), + (np.uint16, np.double), + (np.uint16, np.uint16)) diff --git a/src/silx/math/test/test_histogramnd_nominal.py b/src/silx/math/test/test_histogramnd_nominal.py new file mode 100644 index 0000000..031a772 --- /dev/null +++ b/src/silx/math/test/test_histogramnd_nominal.py @@ -0,0 +1,937 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2021 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +""" +Nominal tests of the histogramnd function. +""" + +import unittest +import pytest + +import numpy as np + +from silx.math.chistogramnd import chistogramnd as histogramnd +from silx.math import Histogramnd + + +def _get_bin_edges(histo_range, n_bins, n_dims): + edges = [] + for i_dim in range(n_dims): + edges.append(histo_range[i_dim, 0] + + np.arange(n_bins[i_dim] + 1) * + (histo_range[i_dim, 1] - histo_range[i_dim, 0]) / + n_bins[i_dim]) + return tuple(edges) + + +# ============================================================== +# ============================================================== +# ============================================================== + + +class _Test_chistogramnd_nominal(unittest.TestCase): + """ + Unit tests of the histogramnd function. + """ + __test__ = False # ignore abstract classe + + ndims = None + + def setUp(self): + if type(self).__name__.startswith("_"): + self.skipTest("Abstract class") + ndims = self.ndims + self.tested_dim = ndims-1 + + if ndims is None: + raise ValueError('ndims class member not set.') + + sample = np.array([5.5, -3.3, + 0., -0.5, + 3.3, 8.8, + -7.7, 6.0, + -4.0]) + + weights = np.array([500.5, -300.3, + 0.01, -0.5, + 300.3, 800.8, + -700.7, 600.6, + -400.4]) + + n_elems = len(sample) + + if ndims == 1: + shape = (n_elems,) + else: + shape = (n_elems, ndims) + + self.sample = np.zeros(shape=shape, dtype=sample.dtype) + if ndims == 1: + self.sample = sample + else: + self.sample[..., ndims-1] = sample + + self.weights = weights + + # the tests are performed along one dimension, + # all the other bins indices along the other dimensions + # are expected to be 2 + # (e.g : when testing a 2D sample : [0, x] will go into + # bin [2, y] because of the bin ranges [-2, 2] and n_bins = 4 + # for the first dimension) + self.other_axes_index = 2 + self.histo_range = np.repeat([[-2., 2.]], ndims, axis=0) + self.histo_range[ndims-1] = [-4., 6.] + + self.n_bins = np.array([4]*ndims) + self.n_bins[ndims-1] = 5 + + if ndims == 1: + def fill_histo(h, v, dim, op=None): + if op: + h[:] = op(h[:], v) + else: + h[:] = v + self.fill_histo = fill_histo + else: + def fill_histo(h, v, dim, op=None): + idx = [self.other_axes_index]*len(h.shape) + idx[dim] = slice(0, None) + idx = tuple(idx) + if op: + h[idx] = op(h[idx], v) + else: + h[idx] = v + self.fill_histo = fill_histo + + def test_nominal(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul, bin_edges = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights) + + expected_edges = _get_bin_edges(self.histo_range, + self.n_bins, + self.ndims) + + self.assertEqual(cumul.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + for i_edges, edges in enumerate(expected_edges): + self.assertTrue(np.array_equal(bin_edges[i_edges], + expected_edges[i_edges]), + msg='Testing bin_edges for dim {0}' + ''.format(i_edges+1)) + + def test_nominal_wh_dtype(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.float32) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul, bin_edges = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + wh_dtype=np.float32) + + self.assertEqual(cumul.dtype, np.float32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.allclose(cumul, expected_c)) + + def test_nominal_uncontiguous_sample(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + shape = list(self.sample.shape) + shape[0] *= 2 + sample = np.zeros(shape, dtype=self.sample.dtype) + uncontig_sample = sample[::2, ...] + uncontig_sample[:] = self.sample + + self.assertFalse(uncontig_sample.flags['C_CONTIGUOUS'], + msg='Making sure the array is not contiguous.') + + histo, cumul, bin_edges = histogramnd(uncontig_sample, + self.histo_range, + self.n_bins, + weights=self.weights) + + self.assertEqual(cumul.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_nominal_uncontiguous_weights(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + shape = list(self.weights.shape) + shape[0] *= 2 + weights = np.zeros(shape, dtype=self.weights.dtype) + uncontig_weights = weights[::2, ...] + uncontig_weights[:] = self.weights + + self.assertFalse(uncontig_weights.flags['C_CONTIGUOUS'], + msg='Making sure the array is not contiguous.') + + histo, cumul, bin_edges = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=uncontig_weights) + + self.assertEqual(cumul.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_nominal_wo_weights(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=None)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(cumul is None) + + def test_nominal_wo_weights_w_cumul(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + + # creating an array of ones just to make sure that + # it is not cleared by histogramnd + cumul_in = np.ones(self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=None, + weighted_histo=cumul_in)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(cumul is None) + self.assertTrue(np.array_equal(cumul_in, + np.ones(shape=self.n_bins, + dtype=np.double))) + + def test_nominal_wo_weights_w_histo(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + + # creating an array of ones just to make sure that + # it is not cleared by histogramnd + histo_in = np.ones(self.n_bins, dtype=np.uint32) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=None, + histo=histo_in)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h + 1)) + self.assertTrue(cumul is None) + self.assertEqual(id(histo), id(histo_in)) + + def test_nominal_last_bin_closed(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 2]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 1101.1]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_int32_weights_double_weights_range(self): + """ + """ + weight_min = -299.9 # ===> will be cast to -299 + weight_max = 499.9 # ===> will be cast to 499 + + expected_h_tpl = np.array([0, 1, 1, 1, 0]) + expected_c_tpl = np.array([0., 0., 0., 300., 0.]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights.astype(np.int32), + weight_min=weight_min, + weight_max=weight_max)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_reuse_histo(self): + """ + """ + + expected_h_tpl = np.array([2, 3, 2, 2, 2]) + expected_c_tpl = np.array([0.0, -7007, -5.0, 0.1, 3003.0]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights)[0:2] + + sample_2 = self.sample[:] + if len(sample_2.shape) == 1: + idx = (slice(0, None),) + else: + idx = slice(0, None), self.tested_dim + + sample_2[idx] += 2 + + histo_2, cumul = histogramnd(sample_2, # <==== !! + self.histo_range, + self.n_bins, + weights=10 * self.weights, # <==== !! + histo=histo)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + self.assertEqual(id(histo), id(histo_2)) + + def test_reuse_cumul(self): + """ + """ + + expected_h_tpl = np.array([0, 2, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -7007.5, -4.99, 300.4, 3503.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights)[0:2] + + sample_2 = self.sample[:] + if len(sample_2.shape) == 1: + idx = (slice(0, None),) + else: + idx = slice(0, None), self.tested_dim + + sample_2[idx] += 2 + + histo, cumul_2 = histogramnd(sample_2, # <==== !! + self.histo_range, + self.n_bins, + weights=10 * self.weights, # <==== !! + weighted_histo=cumul)[0:2] + + self.assertEqual(cumul.dtype, np.float64) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.allclose(cumul, expected_c, rtol=10e-15)) + self.assertEqual(id(cumul), id(cumul_2)) + + def test_reuse_cumul_float(self): + """ + """ + + expected_h_tpl = np.array([0, 2, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -7007.5, -4.99, 300.4, 3503.5], + dtype=np.float32) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights)[0:2] + + # converting the cumul array to float + cumul = cumul.astype(np.float32) + + sample_2 = self.sample[:] + if len(sample_2.shape) == 1: + idx = (slice(0, None),) + else: + idx = slice(0, None), self.tested_dim + + sample_2[idx] += 2 + + histo, cumul_2 = histogramnd(sample_2, # <==== !! + self.histo_range, + self.n_bins, + weights=10 * self.weights, # <==== !! + weighted_histo=cumul)[0:2] + + self.assertEqual(cumul.dtype, np.float32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertEqual(id(cumul), id(cumul_2)) + self.assertTrue(np.allclose(cumul, expected_c, rtol=10e-15)) + +class _Test_Histogramnd_nominal(unittest.TestCase): + """ + Unit tests of the Histogramnd class. + """ + __test__ = False # ignore abstract class + + ndims = None + + def setUp(self): + ndims = self.ndims + if ndims is None: + self.skipTest("Abstract class") + self.tested_dim = ndims-1 + + if ndims is None: + raise ValueError('ndims class member not set.') + + sample = np.array([5.5, -3.3, + 0., -0.5, + 3.3, 8.8, + -7.7, 6.0, + -4.0]) + + weights = np.array([500.5, -300.3, + 0.01, -0.5, + 300.3, 800.8, + -700.7, 600.6, + -400.4]) + + n_elems = len(sample) + + if ndims == 1: + shape = (n_elems,) + else: + shape = (n_elems, ndims) + + self.sample = np.zeros(shape=shape, dtype=sample.dtype) + if ndims == 1: + self.sample = sample + else: + self.sample[..., ndims-1] = sample + + self.weights = weights + + # the tests are performed along one dimension, + # all the other bins indices along the other dimensions + # are expected to be 2 + # (e.g : when testing a 2D sample : [0, x] will go into + # bin [2, y] because of the bin ranges [-2, 2] and n_bins = 4 + # for the first dimension) + self.other_axes_index = 2 + self.histo_range = np.repeat([[-2., 2.]], ndims, axis=0) + self.histo_range[ndims-1] = [-4., 6.] + + self.n_bins = np.array([4]*ndims) + self.n_bins[ndims-1] = 5 + + if ndims == 1: + def fill_histo(h, v, dim, op=None): + if op: + h[:] = op(h[:], v) + else: + h[:] = v + self.fill_histo = fill_histo + else: + def fill_histo(h, v, dim, op=None): + idx = [self.other_axes_index]*len(h.shape) + idx[dim] = slice(0, None) + idx = tuple(idx) + if op: + h[idx] = op(h[idx], v) + else: + h[idx] = v + self.fill_histo = fill_histo + + def test_nominal(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights) + + histo, cumul, bin_edges = histo + + expected_edges = _get_bin_edges(self.histo_range, + self.n_bins, + self.ndims) + + self.assertEqual(cumul.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + for i_edges, edges in enumerate(expected_edges): + self.assertTrue(np.array_equal(bin_edges[i_edges], + expected_edges[i_edges]), + msg='Testing bin_edges for dim {0}' + ''.format(i_edges+1)) + + def test_nominal_wh_dtype(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.float32) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul, bin_edges = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + wh_dtype=np.float32) + + self.assertEqual(cumul.dtype, np.float32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.allclose(cumul, expected_c)) + + def test_nominal_uncontiguous_sample(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + shape = list(self.sample.shape) + shape[0] *= 2 + sample = np.zeros(shape, dtype=self.sample.dtype) + uncontig_sample = sample[::2, ...] + uncontig_sample[:] = self.sample + + self.assertFalse(uncontig_sample.flags['C_CONTIGUOUS'], + msg='Making sure the array is not contiguous.') + + histo, cumul, bin_edges = Histogramnd(uncontig_sample, + self.histo_range, + self.n_bins, + weights=self.weights) + + self.assertEqual(cumul.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_nominal_uncontiguous_weights(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + shape = list(self.weights.shape) + shape[0] *= 2 + weights = np.zeros(shape, dtype=self.weights.dtype) + uncontig_weights = weights[::2, ...] + uncontig_weights[:] = self.weights + + self.assertFalse(uncontig_weights.flags['C_CONTIGUOUS'], + msg='Making sure the array is not contiguous.') + + histo, cumul, bin_edges = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=uncontig_weights) + + self.assertEqual(cumul.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_nominal_wo_weights(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + + histo, cumul = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=None)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(cumul is None) + + def test_nominal_last_bin_closed(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 2]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 1101.1]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_int32_weights_double_weights_range(self): + """ + """ + weight_min = -299.9 # ===> will be cast to -299 + weight_max = 499.9 # ===> will be cast to 499 + + expected_h_tpl = np.array([0, 1, 1, 1, 0]) + expected_c_tpl = np.array([0., 0., 0., 300., 0.]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo, cumul = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights.astype(np.int32), + weight_min=weight_min, + weight_max=weight_max)[0:2] + + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def test_nominal_no_sample(self): + """ + """ + + histo_inst = Histogramnd(None, + self.histo_range, + self.n_bins) + + histo, weighted_histo, edges = histo_inst + + self.assertIsNone(histo) + self.assertIsNone(weighted_histo) + self.assertIsNone(edges) + self.assertIsNone(histo_inst.histo) + self.assertIsNone(histo_inst.weighted_histo) + self.assertIsNone(histo_inst.edges) + + def test_empty_init_accumulate(self): + """ + """ + expected_h_tpl = np.array([2, 1, 1, 1, 1]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo_inst = Histogramnd(None, + self.histo_range, + self.n_bins) + + histo_inst.accumulate(self.sample, + weights=self.weights) + + histo = histo_inst.histo + cumul = histo_inst.weighted_histo + bin_edges = histo_inst.edges + + expected_edges = _get_bin_edges(self.histo_range, + self.n_bins, + self.ndims) + + self.assertEqual(cumul.dtype, np.float64) + self.assertEqual(histo.dtype, np.uint32) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + for i_edges, edges in enumerate(expected_edges): + self.assertTrue(np.array_equal(bin_edges[i_edges], + expected_edges[i_edges]), + msg='Testing bin_edges for dim {0}' + ''.format(i_edges+1)) + + def test_accumulate(self): + """ + """ + + expected_h_tpl = np.array([2, 3, 2, 2, 2]) + expected_c_tpl = np.array([-700.7, -7007.5, -4.99, 300.4, 3503.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo_inst = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights) + + sample_2 = self.sample[:] + if len(sample_2.shape) == 1: + idx = (slice(0, None),) + else: + idx = slice(0, None), self.tested_dim + + sample_2[idx] += 2 + + histo_inst.accumulate(sample_2, # <==== !! + weights=10 * self.weights) # <==== !! + + histo = histo_inst.histo + cumul = histo_inst.weighted_histo + bin_edges = histo_inst.edges + + self.assertEqual(cumul.dtype, np.float64) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.allclose(cumul, expected_c, rtol=10e-15)) + + def test_accumulate_no_weights(self): + """ + """ + + expected_h_tpl = np.array([2, 3, 2, 2, 2]) + expected_c_tpl = np.array([-700.7, -0.5, 0.01, 300.3, 500.5]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo_inst = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights) + + sample_2 = self.sample[:] + if len(sample_2.shape) == 1: + idx = (slice(0, None),) + else: + idx = slice(0, None), self.tested_dim + + sample_2[idx] += 2 + + histo_inst.accumulate(sample_2) # <==== !! + + histo = histo_inst.histo + cumul = histo_inst.weighted_histo + bin_edges = histo_inst.edges + + self.assertEqual(cumul.dtype, np.float64) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.allclose(cumul, expected_c, rtol=10e-15)) + + def test_accumulate_no_weights_at_init(self): + """ + """ + + expected_h_tpl = np.array([2, 3, 2, 2, 2]) + expected_c_tpl = np.array([0.0, -700.7, -0.5, 0.01, 300.3]) + + expected_h = np.zeros(shape=self.n_bins, dtype=np.double) + expected_c = np.zeros(shape=self.n_bins, dtype=np.double) + + self.fill_histo(expected_h, expected_h_tpl, self.ndims-1) + self.fill_histo(expected_c, expected_c_tpl, self.ndims-1) + + histo_inst = Histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=None) # <==== !! + + cumul = histo_inst.weighted_histo + self.assertIsNone(cumul) + + sample_2 = self.sample[:] + if len(sample_2.shape) == 1: + idx = (slice(0, None),) + else: + idx = slice(0, None), self.tested_dim + + sample_2[idx] += 2 + + histo_inst.accumulate(sample_2, + weights=self.weights) # <==== !! + + histo = histo_inst.histo + cumul = histo_inst.weighted_histo + bin_edges = histo_inst.edges + + self.assertEqual(cumul.dtype, np.float64) + self.assertTrue(np.array_equal(histo, expected_h)) + self.assertTrue(np.array_equal(cumul, expected_c)) + + def testNoneNativeTypes(self): + type = self.sample.dtype.newbyteorder("B") + sampleB = self.sample.astype(type) + + type = self.sample.dtype.newbyteorder("L") + sampleL = self.sample.astype(type) + + histo_inst = Histogramnd(sampleB, + self.histo_range, + self.n_bins, + weights=self.weights) + + histo_inst = Histogramnd(sampleL, + self.histo_range, + self.n_bins, + weights=self.weights) + + +class Test_chistogram_nominal_1d(_Test_chistogramnd_nominal): + __test__ = True # because _Test_chistogramnd_nominal is ignored + ndims = 1 + + +class Test_chistogram_nominal_2d(_Test_chistogramnd_nominal): + __test__ = True # because _Test_chistogramnd_nominal is ignored + ndims = 2 + + +class Test_chistogram_nominal_3d(_Test_chistogramnd_nominal): + __test__ = True # because _Test_chistogramnd_nominal is ignored + ndims = 3 + + +class Test_Histogramnd_nominal_1d(_Test_Histogramnd_nominal): + __test__ = True # because _Test_chistogramnd_nominal is ignored + ndims = 1 + + +class Test_Histogramnd_nominal_2d(_Test_Histogramnd_nominal): + __test__ = True # because _Test_chistogramnd_nominal is ignored + ndims = 2 + + +class Test_Histogramnd_nominal_3d(_Test_Histogramnd_nominal): + __test__ = True # because _Test_chistogramnd_nominal is ignored + ndims = 3 diff --git a/src/silx/math/test/test_histogramnd_vs_np.py b/src/silx/math/test/test_histogramnd_vs_np.py new file mode 100644 index 0000000..d6a8d19 --- /dev/null +++ b/src/silx/math/test/test_histogramnd_vs_np.py @@ -0,0 +1,826 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016-2021 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +""" +Tests for the histogramnd function. +Results are compared to numpy's histogramdd. +""" + +import unittest +import operator + +import numpy as np + +from silx.math.chistogramnd import chistogramnd as histogramnd + +# ============================================================== +# ============================================================== +# ============================================================== + +_RTOL_DICT = {np.float64: 10**-13, + np.float32: 10**-5} + +# ============================================================== +# ============================================================== +# ============================================================== + + +def _add_values_to_array_if_missing(array, values, n_values): + max_in_col = np.any(array[:, ...] == values, axis=0) + + if len(array.shape) == 1: + if not max_in_col: + rnd_idx = np.random.randint(0, + high=len(array)-1, + size=(n_values,)) + array[rnd_idx] = values + else: + for i in range(len(max_in_col)): + if not max_in_col[i]: + rnd_idx = np.random.randint(0, + high=len(array)-1, + size=(n_values,)) + array[rnd_idx, i] = values[i] + + +def _get_values_index(array, values, op=operator.lt): + idx = op(array[:, ...], values) + if array.ndim > 1: + idx = np.all(idx, axis=1) + return np.where(idx)[0] + + +def _get_in_range_indices(array, + minvalues, + maxvalues, + minop=operator.ge, + maxop=operator.lt): + idx = np.logical_and(minop(array, minvalues), + maxop(array, maxvalues)) + if array.ndim > 1: + idx = np.all(idx, axis=1) + return np.where(idx)[0] + + +class _TestHistogramnd(unittest.TestCase): + """ + Unit tests of the histogramnd function. + """ + __test__ = False # ignore abstract class + + sample_rng = None + weights_rng = None + n_dims = None + + filter_min = None + filter_max = None + + histo_range = None + n_bins = None + + dtype_sample = None + dtype_weights = None + + def generate_data(self): + + self.longMessage = True + + int_min = 0 + int_max = 100000 + n_elements = 10**5 + + if self.n_dims == 1: + shape = (n_elements,) + else: + shape = (n_elements, self.n_dims,) + + self.rng_state = np.random.get_state() + + self.state_msg = ('Current RNG state :\n' + '{0}'.format(self.rng_state)) + + sample = np.random.randint(int_min, + high=int_max, + size=shape) + + sample = sample.astype(self.dtype_sample) + sample = (self.sample_rng[0] + + (sample-int_min) * + (self.sample_rng[1]-self.sample_rng[0]) / + (int_max-int_min)).astype(self.dtype_sample) + + weights = np.random.randint(int_min, + high=int_max, + size=(n_elements,)) + weights = weights.astype(self.dtype_weights) + weights = (self.weights_rng[0] + + (weights-int_min) * + (self.weights_rng[1]-self.weights_rng[0]) / + (int_max-int_min)).astype(self.dtype_weights) + + # !!!!!!!!!!!!!!!!!!!!!!!!!!!! + # !!!!!!!!!!!!!!!!!!!!!!!!!!!! + # the bins range are cast to the same type as the sample + # in order to get the same results as numpy + # (which doesnt cast the range) + self.histo_range = np.array(self.histo_range).astype(self.dtype_sample) + + # adding some values that are equal to the max + # in order to test the opened/closed last bin + bins_max = [b[1] for b in self.histo_range] + _add_values_to_array_if_missing(sample, + bins_max, + 100) + + # adding some values that are equal to the min weight value + # in order to test the filters + _add_values_to_array_if_missing(weights, + self.weights_rng[0], + 100) + + # adding some values that are equal to the max weight value + # in order to test the filters + _add_values_to_array_if_missing(weights, + self.weights_rng[1], + 100) + + return sample, weights + + def setUp(self): + if type(self).__name__.startswith("_"): + self.skipTest("Abstract class") + self.sample, self.weights = self.generate_data() + self.rtol = _RTOL_DICT.get(self.dtype_weights, None) + + def array_compare(self, ar_a, ar_b): + if self.rtol is None: + return np.array_equal(ar_a, ar_b) + return np.allclose(ar_a, ar_b, self.rtol) + + def test_bin_ranges(self): + """ + + """ + result_c = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True) + + result_np = np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range) + + for i_edges, edges in enumerate(result_c[2]): + # allclose for now until I can try with the latest version (TBD) + # of numpy + self.assertTrue(np.allclose(edges, + result_np[1][i_edges]), + msg='{0}. Testing bin_edges for dim {1}.' + ''.format(self.state_msg, i_edges+1)) + + def test_last_bin_closed(self): + """ + + """ + result_c = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True) + + result_np = np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range) + + result_np_w = np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range, + weights=self.weights) + + # comparing "hits" + hits_cmp = np.array_equal(result_c[0], + result_np[0]) + # comparing weights + weights_cmp = np.array_equal(result_c[1], + result_np_w[0]) + + self.assertTrue(hits_cmp, msg=self.state_msg) + self.assertTrue(weights_cmp, msg=self.state_msg) + + bins_min = [rng[0] for rng in self.histo_range] + bins_max = [rng[1] for rng in self.histo_range] + inrange_idx = _get_in_range_indices(self.sample, + bins_min, + bins_max, + minop=operator.ge, + maxop=operator.le) + + self.assertEqual(result_c[0].sum(), inrange_idx.shape[0], + msg=self.state_msg) + + # we have to sum the weights using the same precision as the + # histogramnd function + weights_sum = self.weights[inrange_idx].astype(result_c[1].dtype).sum() + self.assertTrue(self.array_compare(result_c[1].sum(), weights_sum), + msg=self.state_msg) + + def test_last_bin_open(self): + """ + + """ + result_c = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=False) + + bins_max = [rng[1] for rng in self.histo_range] + filtered_idx = _get_values_index(self.sample, bins_max) + + result_np = np.histogramdd(self.sample[filtered_idx], + bins=self.n_bins, + range=self.histo_range) + + result_np_w = np.histogramdd(self.sample[filtered_idx], + bins=self.n_bins, + range=self.histo_range, + weights=self.weights[filtered_idx]) + + # comparing "hits" + hits_cmp = np.array_equal(result_c[0], result_np[0]) + # comparing weights + weights_cmp = np.array_equal(result_c[1], + result_np_w[0]) + + self.assertTrue(hits_cmp, msg=self.state_msg) + self.assertTrue(weights_cmp, msg=self.state_msg) + + bins_min = [rng[0] for rng in self.histo_range] + bins_max = [rng[1] for rng in self.histo_range] + inrange_idx = _get_in_range_indices(self.sample, + bins_min, + bins_max, + minop=operator.ge, + maxop=operator.lt) + + self.assertEqual(result_c[0].sum(), len(inrange_idx), + msg=self.state_msg) + # we have to sum the weights using the same precision as the + # histogramnd function + weights_sum = self.weights[inrange_idx].astype(result_c[1].dtype).sum() + self.assertTrue(self.array_compare(result_c[1].sum(), weights_sum), + msg=self.state_msg) + + def test_filter_min(self): + """ + + """ + result_c = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True, + weight_min=self.filter_min) + + # !!!!!!!!!!!!!!!!!!!!!!!!!!!! + filter_min = self.dtype_weights(self.filter_min) + + weight_idx = _get_values_index(self.weights, + filter_min, # <------ !!! + operator.ge) + + result_np = np.histogramdd(self.sample[weight_idx], + bins=self.n_bins, + range=self.histo_range) + + result_np_w = np.histogramdd(self.sample[weight_idx], + bins=self.n_bins, + range=self.histo_range, + weights=self.weights[weight_idx]) + + # comparing "hits" + hits_cmp = np.array_equal(result_c[0], + result_np[0]) + # comparing weights + weights_cmp = np.array_equal(result_c[1], result_np_w[0]) + + self.assertTrue(hits_cmp, msg=self.state_msg) + self.assertTrue(weights_cmp, msg=self.state_msg) + + bins_min = [rng[0] for rng in self.histo_range] + bins_max = [rng[1] for rng in self.histo_range] + inrange_idx = _get_in_range_indices(self.sample[weight_idx], + bins_min, + bins_max, + minop=operator.ge, + maxop=operator.le) + + inrange_idx = weight_idx[inrange_idx] + + self.assertEqual(result_c[0].sum(), len(inrange_idx), + msg=self.state_msg) + + # we have to sum the weights using the same precision as the + # histogramnd function + weights_sum = self.weights[inrange_idx].astype(result_c[1].dtype).sum() + self.assertTrue(self.array_compare(result_c[1].sum(), weights_sum), + msg=self.state_msg) + + def test_filter_max(self): + """ + + """ + result_c = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True, + weight_max=self.filter_max) + + # !!!!!!!!!!!!!!!!!!!!!!!!!!!! + filter_max = self.dtype_weights(self.filter_max) + + weight_idx = _get_values_index(self.weights, + filter_max, # <------ !!! + operator.le) + + result_np = np.histogramdd(self.sample[weight_idx], + bins=self.n_bins, + range=self.histo_range) + + result_np_w = np.histogramdd(self.sample[weight_idx], + bins=self.n_bins, + range=self.histo_range, + weights=self.weights[weight_idx]) + + # comparing "hits" + hits_cmp = np.array_equal(result_c[0], + result_np[0]) + # comparing weights + weights_cmp = np.array_equal(result_c[1], result_np_w[0]) + + self.assertTrue(hits_cmp, msg=self.state_msg) + self.assertTrue(weights_cmp, msg=self.state_msg) + + bins_min = [rng[0] for rng in self.histo_range] + bins_max = [rng[1] for rng in self.histo_range] + inrange_idx = _get_in_range_indices(self.sample[weight_idx], + bins_min, + bins_max, + minop=operator.ge, + maxop=operator.le) + + inrange_idx = weight_idx[inrange_idx] + + self.assertEqual(result_c[0].sum(), len(inrange_idx), + msg=self.state_msg) + + # we have to sum the weights using the same precision as the + # histogramnd function + weights_sum = self.weights[inrange_idx].astype(result_c[1].dtype).sum() + self.assertTrue(self.array_compare(result_c[1].sum(), weights_sum), + msg=self.state_msg) + + def test_filter_minmax(self): + """ + + """ + result_c = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True, + weight_min=self.filter_min, + weight_max=self.filter_max) + + # !!!!!!!!!!!!!!!!!!!!!!!!!!!! + filter_min = self.dtype_weights(self.filter_min) + filter_max = self.dtype_weights(self.filter_max) + + weight_idx = _get_in_range_indices(self.weights, + filter_min, # <------ !!! + filter_max, # <------ !!! + minop=operator.ge, + maxop=operator.le) + + result_np = np.histogramdd(self.sample[weight_idx], + bins=self.n_bins, + range=self.histo_range) + + result_np_w = np.histogramdd(self.sample[weight_idx], + bins=self.n_bins, + range=self.histo_range, + weights=self.weights[weight_idx]) + + # comparing "hits" + hits_cmp = np.array_equal(result_c[0], + result_np[0]) + # comparing weights + weights_cmp = np.array_equal(result_c[1], result_np_w[0]) + + self.assertTrue(hits_cmp) + self.assertTrue(weights_cmp) + + bins_min = [rng[0] for rng in self.histo_range] + bins_max = [rng[1] for rng in self.histo_range] + inrange_idx = _get_in_range_indices(self.sample[weight_idx], + bins_min, + bins_max, + minop=operator.ge, + maxop=operator.le) + + inrange_idx = weight_idx[inrange_idx] + + self.assertEqual(result_c[0].sum(), len(inrange_idx), + msg=self.state_msg) + + # we have to sum the weights using the same precision as the + # histogramnd function + weights_sum = self.weights[inrange_idx].astype(result_c[1].dtype).sum() + self.assertTrue(self.array_compare(result_c[1].sum(), weights_sum), + msg=self.state_msg) + + def test_reuse_histo(self): + """ + + """ + result_c_1 = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True) + + result_np_1 = np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range) + + np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range, + weights=self.weights) + + sample_2, weights_2 = self.generate_data() + + result_c_2 = histogramnd(sample_2, + self.histo_range, + self.n_bins, + weights=weights_2, + last_bin_closed=True, + histo=result_c_1[0]) + + result_np_2 = np.histogramdd(sample_2, + bins=self.n_bins, + range=self.histo_range) + + result_np_w_2 = np.histogramdd(sample_2, + bins=self.n_bins, + range=self.histo_range, + weights=weights_2) + + # comparing "hits" + hits_cmp = np.array_equal(result_c_2[0], + result_np_1[0] + + result_np_2[0]) + # comparing weights + weights_cmp = np.array_equal(result_c_2[1], + result_np_w_2[0]) + + self.assertTrue(hits_cmp, msg=self.state_msg) + self.assertTrue(weights_cmp, msg=self.state_msg) + + def test_reuse_cumul(self): + """ + + """ + result_c = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True) + + np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range) + + result_np_w = np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range, + weights=self.weights) + + sample_2, weights_2 = self.generate_data() + + result_c_2 = histogramnd(sample_2, + self.histo_range, + self.n_bins, + weights=weights_2, + last_bin_closed=True, + weighted_histo=result_c[1]) + + result_np_2 = np.histogramdd(sample_2, + bins=self.n_bins, + range=self.histo_range) + + result_np_w_2 = np.histogramdd(sample_2, + bins=self.n_bins, + range=self.histo_range, + weights=weights_2) + + # comparing "hits" + hits_cmp = np.array_equal(result_c_2[0], + result_np_2[0]) + # comparing weights + + self.assertTrue(hits_cmp, msg=self.state_msg) + self.assertTrue(self.array_compare(result_c_2[1], + result_np_w[0] + result_np_w_2[0]), + msg=self.state_msg) + + def test_reuse_cumul_float(self): + """ + + """ + n_bins = np.array(self.n_bins, ndmin=1) + if len(self.sample.shape) == 2: + if len(n_bins) == self.sample.shape[1]: + shp = tuple([x for x in n_bins]) + else: + shp = (self.n_bins,) * self.sample.shape[1] + cumul = np.zeros(shp, dtype=np.float32) + else: + shp = (self.n_bins,) + cumul = np.zeros(shp, dtype=np.float32) + + result_c_1 = histogramnd(self.sample, + self.histo_range, + self.n_bins, + weights=self.weights, + last_bin_closed=True, + weighted_histo=cumul) + + result_np_1 = np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range) + + result_np_w_1 = np.histogramdd(self.sample, + bins=self.n_bins, + range=self.histo_range, + weights=self.weights) + + # comparing "hits" + hits_cmp = np.array_equal(result_c_1[0], + result_np_1[0]) + + self.assertTrue(hits_cmp, msg=self.state_msg) + self.assertEqual(result_c_1[1].dtype, np.float32, msg=self.state_msg) + + bins_min = [rng[0] for rng in self.histo_range] + bins_max = [rng[1] for rng in self.histo_range] + inrange_idx = _get_in_range_indices(self.sample, + bins_min, + bins_max, + minop=operator.ge, + maxop=operator.le) + weights_sum = \ + self.weights[inrange_idx].astype(np.float32).sum(dtype=np.float64) + self.assertTrue(np.allclose(result_c_1[1].sum(dtype=np.float64), + weights_sum), msg=self.state_msg) + self.assertTrue(np.allclose(result_c_1[1].sum(dtype=np.float64), + result_np_w_1[0].sum(dtype=np.float64)), + msg=self.state_msg) + + +class _TestHistogramnd_1d(_TestHistogramnd): + """ + Unit tests of the 1D histogramnd function. + """ + sample_rng = [-55., 100.] + weights_rng = [-70., 150.] + n_dims = 1 + filter_min = -15.6 + filter_max = 85.7 + + histo_range = [[-30.2, 90.3]] + n_bins = 30 + + dtype = None + + +class _TestHistogramnd_2d(_TestHistogramnd): + """ + Unit tests of the 1D histogramnd function. + """ + sample_rng = [-50.2, 100.99] + weights_rng = [70., 150.] + n_dims = 2 + filter_min = 81.7 + filter_max = 135.3 + + histo_range = [[10., 90.], [20., 70.]] + n_bins = 30 + + dtype = None + + +class _TestHistogramnd_3d(_TestHistogramnd): + """ + Unit tests of the 1D histogramnd function. + """ + sample_rng = [10.2, 200.9] + weights_rng = [0., 100.] + n_dims = 3 + filter_min = 31.5 + filter_max = 83.7 + + histo_range = [[30.8, 150.2], [20.1, 90.9], [10.1, 195.]] + n_bins = 30 + + dtype = None + + +# ################################################################ +# ################################################################ +# ################################################################ +# ################################################################ + + +class TestHistogramnd_1d_double_double(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.double + + +class TestHistogramnd_1d_double_float(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.float32 + + +class TestHistogramnd_1d_double_int32(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.int32 + + +class TestHistogramnd_1d_float_double(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.double + + +class TestHistogramnd_1d_float_float(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.float32 + + +class TestHistogramnd_1d_float_int32(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.int32 + + +class TestHistogramnd_1d_int32_double(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.double + + +class TestHistogramnd_1d_int32_float(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.float32 + + +class TestHistogramnd_1d_int32_int32(_TestHistogramnd_1d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.int32 + + +class TestHistogramnd_2d_double_double(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.double + + +class TestHistogramnd_2d_double_float(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.float32 + + +class TestHistogramnd_2d_double_int32(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.int32 + + +class TestHistogramnd_2d_float_double(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.double + + +class TestHistogramnd_2d_float_float(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.float32 + + +class TestHistogramnd_2d_float_int32(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.int32 + + +class TestHistogramnd_2d_int32_double(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.double + + +class TestHistogramnd_2d_int32_float(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.float32 + + +class TestHistogramnd_2d_int32_int32(_TestHistogramnd_2d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.int32 + + +class TestHistogramnd_3d_double_double(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.double + + +class TestHistogramnd_3d_double_float(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.float32 + + +class TestHistogramnd_3d_double_int32(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.double + dtype_weights = np.int32 + + +class TestHistogramnd_3d_float_double(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.double + + +class TestHistogramnd_3d_float_float(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.float32 + + +class TestHistogramnd_3d_float_int32(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.float32 + dtype_weights = np.int32 + + +class TestHistogramnd_3d_int32_double(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.double + + +class TestHistogramnd_3d_int32_float(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.float32 + + +class TestHistogramnd_3d_int32_int32(_TestHistogramnd_3d): + __test__ = True # because _TestHistogramnd is ignored + dtype_sample = np.int32 + dtype_weights = np.int32 diff --git a/src/silx/math/test/test_interpolate.py b/src/silx/math/test/test_interpolate.py new file mode 100644 index 0000000..146449d --- /dev/null +++ b/src/silx/math/test/test_interpolate.py @@ -0,0 +1,125 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (c) 2019 European Synchrotron Radiation Facility +# +# 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 for interpolate module""" + +__authors__ = ["T. Vincent"] +__license__ = "MIT" +__date__ = "11/07/2019" + + +import unittest + +import numpy +try: + from scipy.interpolate import interpn +except ImportError: + interpn = None + +from silx.utils.testutils import ParametricTestCase +from silx.math import interpolate + + +@unittest.skipUnless(interpn is not None, "scipy missing") +class TestInterp3d(ParametricTestCase): + """Test silx.math.interpolate.interp3d""" + + @staticmethod + def ref_interp3d(data, points): + """Reference implementation of interp3d based on scipy + + :param numpy.ndarray data: 3D floating dataset + :param numpy.ndarray points: Array of points of shape (N, 3) + """ + return interpn( + [numpy.arange(dim, dtype=data.dtype) for dim in data.shape], + data, + points, + method='linear') + + def test_random_data(self): + """Test interp3d with random data""" + size = 32 + npoints = 10 + + ref_data = numpy.random.random((size, size, size)) + ref_points = numpy.random.random(npoints*3).reshape(npoints, 3) * (size -1) + + for dtype in (numpy.float32, numpy.float64): + data = ref_data.astype(dtype) + points = ref_points.astype(dtype) + ref_result = self.ref_interp3d(data, points) + + for method in (u'linear', u'linear_omp'): + with self.subTest(method=method): + result = interpolate.interp3d(data, points, method=method) + self.assertTrue(numpy.allclose(ref_result, result)) + + def test_notfinite_data(self): + """Test interp3d with NaN and inf""" + data = numpy.ones((3, 3, 3), dtype=numpy.float64) + data[0, 0, 0] = numpy.nan + data[2, 2, 2] = numpy.inf + points = numpy.array([(0.5, 0.5, 0.5), + (1.5, 1.5, 1.5)]) + + for method in (u'linear', u'linear_omp'): + with self.subTest(method=method): + result = interpolate.interp3d( + data, points, method=method) + self.assertTrue(numpy.isnan(result[0])) + self.assertTrue(result[1] == numpy.inf) + + def test_points_outside(self): + """Test interp3d with points outside the volume""" + data = numpy.ones((4, 4, 4), dtype=numpy.float64) + points = numpy.array([(-0.1, -0.1, -0.1), + (3.1, 3.1, 3.1), + (-0.1, 1., 1.), + (1., 1., 3.1)]) + + for method in (u'linear', u'linear_omp'): + for fill_value in (numpy.nan, 0., -1.): + with self.subTest(method=method): + result = interpolate.interp3d( + data, points, method=method, fill_value=fill_value) + if numpy.isnan(fill_value): + self.assertTrue(numpy.all(numpy.isnan(result))) + else: + self.assertTrue(numpy.all(numpy.equal(result, fill_value))) + + def test_integer_points(self): + """Test interp3d with integer points coord""" + data = numpy.arange(4**3, dtype=numpy.float64).reshape(4, 4, 4) + points = numpy.array([(0., 0., 0.), + (0., 0., 1.), + (2., 3., 0.), + (3., 3., 3.)]) + + ref_result = data[tuple(points.T.astype(numpy.int32))] + + for method in (u'linear', u'linear_omp'): + with self.subTest(method=method): + result = interpolate.interp3d(data, points, method=method) + self.assertTrue(numpy.allclose(ref_result, result)) diff --git a/src/silx/math/test/test_marchingcubes.py b/src/silx/math/test/test_marchingcubes.py new file mode 100644 index 0000000..5e2b193 --- /dev/null +++ b/src/silx/math/test/test_marchingcubes.py @@ -0,0 +1,174 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2016 European Synchrotron Radiation Facility +# +# 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. +# +# ############################################################################*/ +"""Tests of the marchingcubes module""" + +from __future__ import division + +__authors__ = ["T. Vincent"] +__license__ = "MIT" +__date__ = "17/01/2018" + +import unittest + +import numpy + +from silx.utils.testutils import ParametricTestCase + +from silx.math import marchingcubes + + +class TestMarchingCubes(ParametricTestCase): + """Tests of marching cubes""" + + def assertAllClose(self, array1, array2, msg=None, + rtol=1e-05, atol=1e-08): + """Assert that the 2 numpy.ndarrays are almost equal. + + :param str msg: Message to provide when assert fails + :param float rtol: Relative tolerance, see :func:`numpy.allclose` + :param float atol: Absolute tolerance, see :func:`numpy.allclose` + """ + if not numpy.allclose(array1, array2, rtol, atol): + raise self.failureException(msg) + + def test_cube(self): + """Unit tests with a single cube""" + + # No isosurface + cube_zero = numpy.zeros((2, 2, 2), dtype=numpy.float32) + + result = marchingcubes.MarchingCubes(cube_zero, 1.) + self.assertEqual(result.shape, cube_zero.shape) + self.assertEqual(result.isolevel, 1.) + self.assertEqual(result.invert_normals, True) + + vertices, normals, indices = result + self.assertEqual(len(vertices), 0) + self.assertEqual(len(normals), 0) + self.assertEqual(len(indices), 0) + + # Cube array dimensions: shape = (dim 0, dim 1, dim2) + # + # dim 0 (Z) + # ^ + # | + # 4 +------+ 5 + # /| /| + # / | / | + # 6 +------+ 7| + # | | | | + # |0 +---|--+ 1 -> dim 2 (X) + # | / | / + # |/ |/ + # 2 +------+ 3 + # / + # dim 1 (Y) + + # isosurface perpendicular to dim 0 (Z) + cube = numpy.array( + (((0., 0.), (0., 0.)), + ((1., 1.), (1., 1.))), dtype=numpy.float32) + level = 0.5 + vertices, normals, indices = marchingcubes.MarchingCubes( + cube, level, invert_normals=False) + self.assertAllClose(vertices[:, 0], level) + self.assertAllClose(normals, (1., 0., 0.)) + self.assertEqual(len(indices), 2) + + # isosurface perpendicular to dim 1 (Y) + cube = numpy.array( + (((0., 0.), (1., 1.)), + ((0., 0.), (1., 1.))), dtype=numpy.float32) + level = 0.2 + vertices, normals, indices = marchingcubes.MarchingCubes(cube, level) + self.assertAllClose(vertices[:, 1], level) + self.assertAllClose(normals, (0., -1., 0.)) + self.assertEqual(len(indices), 2) + + # isosurface perpendicular to dim 2 (X) + cube = numpy.array( + (((0., 1.), (0., 1.)), + ((0., 1.), (0., 1.))), dtype=numpy.float32) + level = 0.9 + vertices, normals, indices = marchingcubes.MarchingCubes( + cube, level, invert_normals=False) + self.assertAllClose(vertices[:, 2], level) + self.assertAllClose(normals, (0., 0., 1.)) + self.assertEqual(len(indices), 2) + + # isosurface normal in dim1, dim 0 (Y, Z) plane + cube = numpy.array( + (((0., 0.), (0., 0.)), + ((0., 0.), (1., 1.))), dtype=numpy.float32) + level = 0.5 + vertices, normals, indices = marchingcubes.MarchingCubes(cube, level) + self.assertAllClose(normals[:, 2], 0.) + self.assertEqual(len(indices), 2) + + def test_sampling(self): + """Test different sampling, comparing to reference without sampling""" + isolevel = 0.5 + size = 9 + chessboard = numpy.zeros((size, size, size), dtype=numpy.float32) + chessboard.reshape(-1)[::2] = 1 # OK as long as dimensions are odd + + ref_result = marchingcubes.MarchingCubes(chessboard, isolevel) + + samplings = [ + (2, 1, 1), + (1, 2, 1), + (1, 1, 2), + (2, 2, 2), + (3, 3, 3), + (1, 3, 1), + (1, 1, 3), + ] + + for sampling in samplings: + with self.subTest(sampling=sampling): + sampling = numpy.array(sampling) + + data = 1e6 * numpy.ones( + sampling * size, dtype=numpy.float32) + # Copy ref chessboard in data according to sampling + data[::sampling[0], ::sampling[1], ::sampling[2]] = chessboard + + result = marchingcubes.MarchingCubes(data, isolevel, + sampling=sampling) + # Compare vertices normalized with shape + self.assertAllClose( + ref_result.get_vertices() / ref_result.shape, + result.get_vertices() / result.shape, + atol=0., rtol=0.) + + # Compare normals + # This comparison only works for normals aligned with axes + # otherwise non uniform sampling would make different normals + self.assertAllClose(ref_result.get_normals(), + result.get_normals(), + atol=0., rtol=0.) + + self.assertAllClose(ref_result.get_indices(), + result.get_indices(), + atol=0., rtol=0.) |