diff --git a/nifty2go/spaces/power_space/power_space.py b/nifty2go/spaces/power_space/power_space.py index 67ba4255136f17e076ce05d3f8fe77e752951b55..315ce1ded7dbfdd9e9d054c1ef111f8097f626ab 100644 --- a/nifty2go/spaces/power_space/power_space.py +++ b/nifty2go/spaces/power_space/power_space.py @@ -113,6 +113,35 @@ class PowerSpace(Space): np.log(float(last_bound)), nbin-1, base=np.e) + @staticmethod + def useful_binbounds(space, logarithmic, nbin=None): + if not (isinstance(space, Space) and space.harmonic): + raise ValueError("first argument must be a harmonic space.") + if logarithmic is None and nbin is None: + return None + nbin = None if nbin is None else int(nbin) + logarithmic = bool(logarithmic) + dists = space.get_unique_distances() + if len(dists) < 3: + raise ValueError("Space does not have enough unique k lengths") + lbound = 0.5*(dists[0]+dists[1]) + rbound = 0.5*(dists[-2]+dists[-1]) + dists[0] = lbound + dists[-1] = rbound + if logarithmic: + dists = np.log(dists) + binsz_min = np.max(np.diff(dists)) + nbin_max = int((dists[-1]-dists[0])/binsz_min)+2 + if nbin is None: + nbin = nbin_max + assert nbin >= 3, "nbin must be at least 3" + if nbin > nbin_max: + raise ValueError("nbin is too large") + if logarithmic: + return PowerSpace.logarithmic_binbounds(nbin, lbound, rbound) + else: + return PowerSpace.linear_binbounds(nbin, lbound, rbound) + def __init__(self, harmonic_partner, binbounds=None): super(PowerSpace, self).__init__() self._ignore_for_hash += ['_pindex', '_kindex', '_rho'] @@ -179,7 +208,8 @@ class PowerSpace(Space): binbounds=self._binbounds) def weight(self): - return 1. + # MR FIXME: this will probably change to 1 soon + return np.asarray(self.rho, dtype=np.float64) def get_distance_array(self): return self.kindex.copy() diff --git a/test/common.py b/test/common.py index 83a9bb6c61ceca990592003b4792c95b7e8bd295..c71bf0b844af614892ccaa9fc08ae9f630bc17ae 100644 --- a/test/common.py +++ b/test/common.py @@ -44,35 +44,3 @@ def generate_spaces(): def generate_harmonic_spaces(): spaces = [RGSpace(4, harmonic=True), LMSpace(5)] return spaces - - -def marco_binbounds(space, logarithmic, nbin=None): - """Only for testing purposes. DO NOT USE IN REAL LIFE!""" - if logarithmic is None and nbin is None: - return None - if not (isinstance(space, Space) and space.harmonic): - raise ValueError("space must be a harmonic space.") - logarithmic = bool(logarithmic) - if nbin is not None: - nbin = int(nbin) - assert nbin >= 3, "nbin must be at least 3" - # equidistant binning (linear or log) - # MR FIXME: this needs to improve - kindex = space.get_unique_distances() - if (logarithmic): - k = np.r_[0, np.log(kindex[1:])] - else: - k = kindex - dk = np.max(k[2:] - k[1:-1]) # minimum dk to avoid empty bins - if(nbin is None): - nbin = int((k[-1] - 0.5 * (k[2] + k[1])) / - dk - 0.5) # maximal nbin - else: - nbin = min(int(nbin), int( - (k[-1] - 0.5 * (k[2] + k[1])) / dk + 2.5)) - dk = (k[-1] - 0.5 * (k[2] + k[1])) / (nbin - 2.5) - bb = np.r_[0.5 * (3 * k[1] - k[2]), - 0.5 * (k[1] + k[2]) + dk * np.arange(nbin-2)] - if(logarithmic): - bb = np.exp(bb) - return bb diff --git a/test/test_operators/test_laplace_operator.py b/test/test_operators/test_laplace_operator.py index 22e7505482a3c9be9679f09a0e5afb00fec49832..b6eea91b772474fb09968f8eaf6371bac4a56067 100644 --- a/test/test_operators/test_laplace_operator.py +++ b/test/test_operators/test_laplace_operator.py @@ -21,14 +21,15 @@ import numpy as np import nifty2go as ift from numpy.testing import assert_allclose from itertools import product -from test.common import expand, marco_binbounds +from test.common import expand class LaplaceOperatorTests(unittest.TestCase): @expand(product([None, False, True], [False, True], [10, 100, 1000])) def test_Laplace(self, log1, log2, sz): s = ift.RGSpace(sz, harmonic=True) - p = ift.PowerSpace(s, binbounds=marco_binbounds(s, logarithmic=log1)) + bb = ift.PowerSpace.useful_binbounds(s, logarithmic=log1) + p = ift.PowerSpace(s, binbounds=bb) L = ift.LaplaceOperator(p, logarithmic=log2) arr = np.random.random(p.shape[0]) fp = ift.Field(p, val=arr) diff --git a/test/test_operators/test_smoothing_operator.py b/test/test_operators/test_smoothing_operator.py index fd05868a95dc4e250bdf853ace611674fd07ebe5..1b908bc0a9ff4df28463149f169808ce9b38ff14 100644 --- a/test/test_operators/test_smoothing_operator.py +++ b/test/test_operators/test_smoothing_operator.py @@ -27,7 +27,7 @@ from nifty2go import Field,\ DirectSmoothingOperator from itertools import product -from test.common import expand, marco_binbounds +from test.common import expand def _get_rtol(tp): @@ -94,8 +94,8 @@ class SmoothingOperator_Tests(unittest.TestCase): def test_smooth_irregular1(self, sz, log, sigma, tp): tol = _get_rtol(tp) sp = RGSpace(sz, harmonic=True) - ps = PowerSpace(sp, binbounds=marco_binbounds( - sp, logarithmic=log, nbin=sz)) + bb = PowerSpace.useful_binbounds(sp, logarithmic=log) + ps = PowerSpace(sp, binbounds=bb) smo = DirectSmoothingOperator(ps, sigma=sigma) inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4, dtype=tp) @@ -107,7 +107,8 @@ class SmoothingOperator_Tests(unittest.TestCase): def test_smooth_irregular2(self, sz1, sz2, log, sigma, tp): tol = _get_rtol(tp) sp = RGSpace([sz1, sz2], harmonic=True) - ps = PowerSpace(sp, binbounds=marco_binbounds(sp, logarithmic=log)) + bb = PowerSpace.useful_binbounds(sp, logarithmic=log) + ps = PowerSpace(sp, binbounds=bb) smo = DirectSmoothingOperator(ps, sigma=sigma) inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4, dtype=tp) diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py index 6e31386320675aa70ba01a3f7e9d60968df01381..68ad9ac3b94c4049e1153186b716fda926a7b819 100644 --- a/test/test_spaces/test_power_space.py +++ b/test/test_spaces/test_power_space.py @@ -24,7 +24,7 @@ import numpy as np from numpy.testing import assert_, assert_equal, assert_almost_equal,\ assert_raises from nifty2go import PowerSpace, RGSpace, Space, LMSpace -from test.common import expand, marco_binbounds +from test.common import expand from itertools import product, chain HARMONIC_SPACES = [RGSpace((8,), harmonic=True), @@ -46,7 +46,7 @@ CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT, # [harmonic_partner, logarithmic, nbin, binbounds, expected] CONSTRUCTOR_CONFIGS = [ - [1, False, None, None, {'error': ValueError}], + [1, False, None, None, {'error': (ValueError, NotImplementedError)}], [RGSpace((8,)), False, None, None, {'error': ValueError}], [RGSpace((8,), harmonic=True), None, None, None, { 'harmonic': True, @@ -61,14 +61,14 @@ CONSTRUCTOR_CONFIGS = [ }], [RGSpace((8,), harmonic=True), True, None, None, { 'harmonic': True, - 'shape': (2,), - 'dim': 2, + 'shape': (4,), + 'dim': 4, 'total_volume': 8.0, 'harmonic_partner': RGSpace((8,), harmonic=True), - 'binbounds': (0.70710678118654757,), - 'pindex': np.array([0, 1, 1, 1, 1, 1, 1, 1]), - 'kindex': np.array([0., 2.28571429]), - 'rho': np.array([1, 7]), + 'binbounds': (0.5, 1.3228756555322954, 3.5), + 'pindex': np.array([0, 1, 2, 2, 3, 2, 2, 1]), + 'kindex': np.array([0., 1., 2.5, 4.]), + 'rho': np.array([1, 2, 4, 1]), }], ] @@ -98,9 +98,8 @@ class PowerSpaceConsistencyCheck(unittest.TestCase): @expand(CONSISTENCY_CONFIGS) def test_rhopindexConsistency(self, harmonic_partner, binbounds, nbin, logarithmic): - p = PowerSpace(harmonic_partner=harmonic_partner, - binbounds=marco_binbounds(harmonic_partner, - logarithmic, nbin)) + bb = PowerSpace.useful_binbounds(harmonic_partner, logarithmic, nbin) + p = PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb) assert_equal(np.bincount(p.pindex.flatten()), p.rho, err_msg='rho is not equal to pindex degeneracy') @@ -112,13 +111,13 @@ class PowerSpaceFunctionalityTest(unittest.TestCase): logarithmic, nbin, binbounds, expected): if 'error' in expected: with assert_raises(expected['error']): - PowerSpace(harmonic_partner=harmonic_partner, - binbounds=marco_binbounds(harmonic_partner, - logarithmic, nbin)) + bb = PowerSpace.useful_binbounds(harmonic_partner, + logarithmic, nbin) + PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb) else: - p = PowerSpace(harmonic_partner=harmonic_partner, - binbounds=marco_binbounds(harmonic_partner, - logarithmic, nbin)) + bb = PowerSpace.useful_binbounds(harmonic_partner, + logarithmic, nbin) + p = PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb) for key, value in expected.items(): if isinstance(value, np.ndarray): assert_almost_equal(getattr(p, key), value) @@ -132,4 +131,4 @@ class PowerSpaceFunctionalityTest(unittest.TestCase): def test_weight(self): p = PowerSpace(harmonic_partner=RGSpace(10,harmonic=True)) - assert_almost_equal(p.weight(),1.) + assert_almost_equal(p.weight(),p.rho)