diff options
Diffstat (limited to 'silx/math/test')
-rw-r--r-- | silx/math/test/__init__.py | 51 | ||||
-rw-r--r-- | silx/math/test/benchmark_combo.py | 204 | ||||
-rw-r--r-- | silx/math/test/histo_benchmarks.py | 269 | ||||
-rw-r--r-- | silx/math/test/test_HistogramndLut_nominal.py | 571 | ||||
-rw-r--r-- | silx/math/test/test_combo.py | 168 | ||||
-rw-r--r-- | silx/math/test/test_histogramnd_error.py | 512 | ||||
-rw-r--r-- | silx/math/test/test_histogramnd_nominal.py | 933 | ||||
-rw-r--r-- | silx/math/test/test_histogramnd_vs_np.py | 848 | ||||
-rw-r--r-- | silx/math/test/test_marchingcubes.py | 188 |
9 files changed, 3744 insertions, 0 deletions
diff --git a/silx/math/test/__init__.py b/silx/math/test/__init__.py new file mode 100644 index 0000000..7171bda --- /dev/null +++ b/silx/math/test/__init__.py @@ -0,0 +1,51 @@ +# 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__ = "04/07/2016" + +import unittest + +from .test_histogramnd_error import suite as test_histo_error +from .test_histogramnd_nominal import suite as test_histo_nominal +from .test_histogramnd_vs_np import suite as test_histo_vs_np +from .test_HistogramndLut_nominal import suite as test_histolut_nominal +from ..fit.test import suite as test_fit_suite +from .test_marchingcubes import suite as test_marchingcubes_suite +from ..medianfilter.test import suite as test_medianfilter_suite +from .test_combo import suite as test_combo_suite + + +def suite(): + test_suite = unittest.TestSuite() + test_suite.addTest(test_histo_nominal()) + test_suite.addTest(test_histo_error()) + test_suite.addTest(test_histo_vs_np()) + test_suite.addTest(test_fit_suite()) + test_suite.addTest(test_histolut_nominal()) + test_suite.addTest(test_marchingcubes_suite()) + test_suite.addTest(test_medianfilter_suite()) + test_suite.addTest(test_combo_suite()) + return test_suite diff --git a/silx/math/test/benchmark_combo.py b/silx/math/test/benchmark_combo.py new file mode 100644 index 0000000..4acc26b --- /dev/null +++ b/silx/math/test/benchmark_combo.py @@ -0,0 +1,204 @@ +# 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__ = "27/03/2017" + + +import logging +import os.path +import time +import unittest + +import numpy + +from silx.test.utils import ParametricTestCase, temp_dir + +from silx.math import combo + + +logging.basicConfig() +_logger = logging.getLogger(__name__) +_logger.setLevel(logging.DEBUG) + + +class BenchmarkMinMax(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() + + +def suite(): + test_suite = unittest.TestSuite() + test_suite.addTests( + unittest.defaultTestLoader.loadTestsFromTestCase(BenchmarkMinMax)) + return test_suite + + +if __name__ == '__main__': + unittest.main(defaultTest="suite") diff --git a/silx/math/test/histo_benchmarks.py b/silx/math/test/histo_benchmarks.py new file mode 100644 index 0000000..7d3216d --- /dev/null +++ b/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/silx/math/test/test_HistogramndLut_nominal.py b/silx/math/test/test_HistogramndLut_nominal.py new file mode 100644 index 0000000..9c356bd --- /dev/null +++ b/silx/math/test/test_HistogramndLut_nominal.py @@ -0,0 +1,571 @@ +# 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. +# +# ############################################################################*/ +""" +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. + """ + + ndims = None + + def setUp(self): + 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) + 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)) + + +class TestHistogramndLut_nominal_1d(_TestHistogramndLut_nominal): + ndims = 1 + + +class TestHistogramndLut_nominal_2d(_TestHistogramndLut_nominal): + ndims = 2 + + +class TestHistogramndLut_nominal_3d(_TestHistogramndLut_nominal): + ndims = 3 + + +# ============================================================== +# ============================================================== +# ============================================================== + + +test_cases = (TestHistogramndLut_nominal_1d, + TestHistogramndLut_nominal_2d, + TestHistogramndLut_nominal_3d,) + + +def suite(): + loader = unittest.defaultTestLoader + test_suite = unittest.TestSuite() + for test_class in test_cases: + tests = loader.loadTestsFromTestCase(test_class) + test_suite.addTests(tests) + return test_suite + +if __name__ == '__main__': + unittest.main(defaultTest="suite") diff --git a/silx/math/test/test_combo.py b/silx/math/test/test_combo.py new file mode 100644 index 0000000..92eb8b4 --- /dev/null +++ b/silx/math/test/test_combo.py @@ -0,0 +1,168 @@ +# 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 combo module""" + +from __future__ import division + +__authors__ = ["T. Vincent"] +__license__ = "MIT" +__date__ = "20/12/2016" + + +import unittest + +import numpy + +from silx.test.utils import ParametricTestCase + +from silx.math.combo import min_max + + +class TestMinMax(ParametricTestCase): + """Tests of min max combo""" + + FLOATING_DTYPES = 'float32', 'float64' + SIGNED_INT_DTYPES = 'uint8', 'uint16', 'uint32', 'uint64' + UNSIGNED_INT_DTYPES = 'uint8', 'uint16', 'uint32', 'uint64' + DTYPES = FLOATING_DTYPES + SIGNED_INT_DTYPES + UNSIGNED_INT_DTYPES + + def _test_min_max(self, data, min_positive): + """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 + """ + result = min_max(data, min_positive) + + minimum = numpy.nanmin(data) + if numpy.isnan(minimum): # All NaNs + self.assertTrue(numpy.isnan(result.minimum)) + self.assertEqual(result.argmin, 0) + + else: + self.assertEqual(result.minimum, minimum) + + argmin = numpy.where(data == minimum)[0][0] + self.assertEqual(result.argmin, argmin) + + maximum = numpy.nanmax(data) + if numpy.isnan(maximum): # All NaNs + self.assertTrue(numpy.isnan(result.maximum)) + self.assertEqual(result.argmax, 0) + + else: + self.assertEqual(result.maximum, maximum) + + argmax = numpy.where(data == maximum)[0][0] + self.assertEqual(result.argmax, argmax) + + if min_positive: + pos_data = data[data > 0] + if len(pos_data) > 0: + min_pos = numpy.min(pos_data) + argmin_pos = numpy.where(data == min_pos)[0][0] + else: + min_pos = None + argmin_pos = None + self.assertEqual(result.min_positive, min_pos) + self.assertEqual(result.argmin_positive, argmin_pos) + + 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) + + def test_nandata(self): + """Test min_max with NaN in data""" + tests = [ + (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 + ] + + 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) + + def test_infdata(self): + """Test min_max with inf.""" + tests = [ + [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 + ] + + 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) + + +def suite(): + test_suite = unittest.TestSuite() + test_suite.addTests( + unittest.defaultTestLoader.loadTestsFromTestCase(TestMinMax)) + return test_suite + + +if __name__ == '__main__': + unittest.main(defaultTest="suite") diff --git a/silx/math/test/test_histogramnd_error.py b/silx/math/test/test_histogramnd_error.py new file mode 100644 index 0000000..575eaf0 --- /dev/null +++ b/silx/math/test/test_histogramnd_error.py @@ -0,0 +1,512 @@ +# 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. + """ + def setUp(self): + raise NotImplementedError('') + + 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. + """ + + 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_errors(_Test_chistogramnd_errors): + """ + Unit tests of the 3D histogramnd error cases. + """ + + 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)) +# ============================================================== +# ============================================================== +# ============================================================== + + +test_cases = (Test_chistogramnd_1D_errors, + Test_chistogramnd_ND_errors,) + + +def suite(): + loader = unittest.defaultTestLoader + test_suite = unittest.TestSuite() + for test_class in test_cases: + tests = loader.loadTestsFromTestCase(test_class) + test_suite.addTests(tests) + return test_suite + +if __name__ == '__main__': + unittest.main(defaultTest="suite") diff --git a/silx/math/test/test_histogramnd_nominal.py b/silx/math/test/test_histogramnd_nominal.py new file mode 100644 index 0000000..f4550e4 --- /dev/null +++ b/silx/math/test/test_histogramnd_nominal.py @@ -0,0 +1,933 @@ +# 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. +# +# ############################################################################*/ +""" +Nominal tests of the histogramnd function. +""" + +import unittest + +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. + """ + + ndims = None + + def setUp(self): + 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) + 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. + """ + + ndims = None + + def setUp(self): + 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) + 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)) + +class Test_chistogram_nominal_1d(_Test_chistogramnd_nominal): + ndims = 1 + + +class Test_chistogram_nominal_2d(_Test_chistogramnd_nominal): + ndims = 2 + + +class Test_chistogram_nominal_3d(_Test_chistogramnd_nominal): + ndims = 3 + + +class Test_Histogramnd_nominal_1d(_Test_Histogramnd_nominal): + ndims = 1 + + +class Test_Histogramnd_nominal_2d(_Test_Histogramnd_nominal): + ndims = 2 + + +class Test_Histogramnd_nominal_3d(_Test_Histogramnd_nominal): + ndims = 3 + + +# ============================================================== +# ============================================================== +# ============================================================== + + +test_cases = (Test_chistogram_nominal_1d, + Test_chistogram_nominal_2d, + Test_chistogram_nominal_3d, + Test_Histogramnd_nominal_1d, + # Test_Histogramnd_nominal_2d, + # Test_Histogramnd_nominal_3d + ) + + +def suite(): + loader = unittest.defaultTestLoader + test_suite = unittest.TestSuite() + for test_class in test_cases: + tests = loader.loadTestsFromTestCase(test_class) + test_suite.addTests(tests) + return test_suite + +if __name__ == '__main__': + unittest.main(defaultTest="suite") diff --git a/silx/math/test/test_histogramnd_vs_np.py b/silx/math/test/test_histogramnd_vs_np.py new file mode 100644 index 0000000..36d71f9 --- /dev/null +++ b/silx/math/test/test_histogramnd_vs_np.py @@ -0,0 +1,848 @@ +# 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 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. + """ + 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): + 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): + dtype_sample = np.double + dtype_weights = np.double + + +class TestHistogramnd_1d_double_float(_TestHistogramnd_1d): + dtype_sample = np.double + dtype_weights = np.float32 + + +class TestHistogramnd_1d_double_int32(_TestHistogramnd_1d): + dtype_sample = np.double + dtype_weights = np.int32 + + +class TestHistogramnd_1d_float_double(_TestHistogramnd_1d): + dtype_sample = np.float32 + dtype_weights = np.double + + +class TestHistogramnd_1d_float_float(_TestHistogramnd_1d): + dtype_sample = np.float32 + dtype_weights = np.float32 + + +class TestHistogramnd_1d_float_int32(_TestHistogramnd_1d): + dtype_sample = np.float32 + dtype_weights = np.int32 + + +class TestHistogramnd_1d_int32_double(_TestHistogramnd_1d): + dtype_sample = np.int32 + dtype_weights = np.double + + +class TestHistogramnd_1d_int32_float(_TestHistogramnd_1d): + dtype_sample = np.int32 + dtype_weights = np.float32 + + +class TestHistogramnd_1d_int32_int32(_TestHistogramnd_1d): + dtype_sample = np.int32 + dtype_weights = np.int32 + + +class TestHistogramnd_2d_double_double(_TestHistogramnd_2d): + dtype_sample = np.double + dtype_weights = np.double + + +class TestHistogramnd_2d_double_float(_TestHistogramnd_2d): + dtype_sample = np.double + dtype_weights = np.float32 + + +class TestHistogramnd_2d_double_int32(_TestHistogramnd_2d): + dtype_sample = np.double + dtype_weights = np.int32 + + +class TestHistogramnd_2d_float_double(_TestHistogramnd_2d): + dtype_sample = np.float32 + dtype_weights = np.double + + +class TestHistogramnd_2d_float_float(_TestHistogramnd_2d): + dtype_sample = np.float32 + dtype_weights = np.float32 + + +class TestHistogramnd_2d_float_int32(_TestHistogramnd_2d): + dtype_sample = np.float32 + dtype_weights = np.int32 + + +class TestHistogramnd_2d_int32_double(_TestHistogramnd_2d): + dtype_sample = np.int32 + dtype_weights = np.double + + +class TestHistogramnd_2d_int32_float(_TestHistogramnd_2d): + dtype_sample = np.int32 + dtype_weights = np.float32 + + +class TestHistogramnd_2d_int32_int32(_TestHistogramnd_2d): + dtype_sample = np.int32 + dtype_weights = np.int32 + + +class TestHistogramnd_3d_double_double(_TestHistogramnd_3d): + dtype_sample = np.double + dtype_weights = np.double + + +class TestHistogramnd_3d_double_float(_TestHistogramnd_3d): + dtype_sample = np.double + dtype_weights = np.float32 + + +class TestHistogramnd_3d_double_int32(_TestHistogramnd_3d): + dtype_sample = np.double + dtype_weights = np.int32 + + +class TestHistogramnd_3d_float_double(_TestHistogramnd_3d): + dtype_sample = np.float32 + dtype_weights = np.double + + +class TestHistogramnd_3d_float_float(_TestHistogramnd_3d): + dtype_sample = np.float32 + dtype_weights = np.float32 + + +class TestHistogramnd_3d_float_int32(_TestHistogramnd_3d): + dtype_sample = np.float32 + dtype_weights = np.int32 + + +class TestHistogramnd_3d_int32_double(_TestHistogramnd_3d): + dtype_sample = np.int32 + dtype_weights = np.double + + +class TestHistogramnd_3d_int32_float(_TestHistogramnd_3d): + dtype_sample = np.int32 + dtype_weights = np.float32 + + +class TestHistogramnd_3d_int32_int32(_TestHistogramnd_3d): + dtype_sample = np.int32 + dtype_weights = np.int32 + + +# ============================================================== +# ============================================================== +# ============================================================== + + +test_cases = (TestHistogramnd_1d_double_double, + TestHistogramnd_1d_double_float, + TestHistogramnd_1d_double_int32, + TestHistogramnd_1d_float_double, + TestHistogramnd_1d_float_float, + TestHistogramnd_1d_float_int32, + TestHistogramnd_1d_int32_double, + TestHistogramnd_1d_int32_float, + TestHistogramnd_1d_int32_int32, + TestHistogramnd_2d_double_double, + TestHistogramnd_2d_double_float, + TestHistogramnd_2d_double_int32, + TestHistogramnd_2d_float_double, + TestHistogramnd_2d_float_float, + TestHistogramnd_2d_float_int32, + TestHistogramnd_2d_int32_double, + TestHistogramnd_2d_int32_float, + TestHistogramnd_2d_int32_int32, + TestHistogramnd_3d_double_double, + TestHistogramnd_3d_double_float, + TestHistogramnd_3d_double_int32, + TestHistogramnd_3d_float_double, + TestHistogramnd_3d_float_float, + TestHistogramnd_3d_float_int32, + TestHistogramnd_3d_int32_double, + TestHistogramnd_3d_int32_float, + TestHistogramnd_3d_int32_int32,) + + +def suite(): + loader = unittest.defaultTestLoader + test_suite = unittest.TestSuite() + for test_class in test_cases: + tests = loader.loadTestsFromTestCase(test_class) + test_suite.addTests(tests) + return test_suite + +if __name__ == '__main__': + unittest.main(defaultTest="suite") diff --git a/silx/math/test/test_marchingcubes.py b/silx/math/test/test_marchingcubes.py new file mode 100644 index 0000000..d6aa8e9 --- /dev/null +++ b/silx/math/test/test_marchingcubes.py @@ -0,0 +1,188 @@ +# 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__ = "05/12/2016" + +import unittest + +import numpy + +from silx.test.utils 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.) + + +test_cases = (TestMarchingCubes,) + + +def suite(): + test_suite = unittest.TestSuite() + for test_class in test_cases: + test_suite.addTests( + unittest.defaultTestLoader.loadTestsFromTestCase(test_class)) + return test_suite + +if __name__ == '__main__': + unittest.main(defaultTest="suite") |