Commit f172caff authored by Martin Reinecke's avatar Martin Reinecke
Browse files

new weight() interface for DomainObjects

parent 18dd49ed
......@@ -127,37 +127,16 @@ class DomainObject(with_metaclass(
"There is no generic dim for DomainObject.")
@abc.abstractmethod
def weight(self, x, power=1, axes=None, inplace=False):
""" Weights the field on this domain with the space's volume-weights.
def weight(self):
""" Returns the volume factors of this domain, either as a floating
point scalar (if the volume factors are all identical) or as a
floating point array with a shape of `self.shape`.
Weights hereby refer to integration weights, as they appear in
discretized integrals. Per default, this function mutliplies each bin
of the field x by its volume, which lets it behave like a density
(top form). However, different powers of the volume can be applied
with the power parameter. The axes parameter specifies which of the
field array's indices correspond to this domain.
Parameters
----------
x : numpy.ndarray
The fields data array.
power : int, *optional*
The power to which the volume-weight is raised (default: 1).
axes : {int, tuple}, *optional*
Specifies the axes of x which represent this domain
(default: None).
If axes==None:
weighting is applied with respect to all axes
inplace : bool, *optional*
If this is True, the weighting is done on the values of x,
if it is False, x is not modified and this method returns a
weighted copy of x (default: False).
Returns
-------
numpy.ndarray
A weighted version of x, with volume-weights raised to the
given power.
float or numpy.ndarray(dtype=float)
Volume factors
Raises
------
......@@ -166,4 +145,4 @@ class DomainObject(with_metaclass(
"""
raise NotImplementedError(
"There is no generic weight-method for DomainObject.")
"There is no generic weight method for DomainObject.")
......@@ -685,20 +685,25 @@ class Field(object):
The weighted field.
"""
new_val = self.get_val(copy=False)
new_val = self.get_val(copy=not inplace)
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
spaces = list(range(len(self.domain)))
for ind, sp in enumerate(self.domain):
if ind in spaces:
new_val = sp.weight(new_val,
power=power,
axes=self.domain_axes[ind],
inplace=inplace)
# we need at most one copy, the rest can happen in place
inplace = True
fct = 1.
for ind in spaces:
wgt = self.domain[ind].weight()
if np.isscalar(wgt):
fct *= wgt
else:
new_shape = np.ones(len(self.shape), dtype=np.int)
new_shape[self.domain_axes[ind][0]:self.domain_axes[ind][-1]+1]=wgt.shape
wgt = wgt.reshape(new_shape)
new_val *= wgt**power
fct = fct**power
if fct!=1:
new_val *= fct
return Field(self.domain, new_val, self.dtype)
......
......@@ -20,5 +20,5 @@ from ..domain_object import DomainObject
class FieldType(DomainObject):
def weight(self, x, power=1, axes=None, inplace=False):
return x if inplace else x.copy()
def weight(self):
return 1.
......@@ -79,8 +79,7 @@ def cast_axis_to_tuple(axis, length=None):
if np.isscalar(axis):
axis = (int(axis),)
else:
raise TypeError(
"Could not convert axis-input to tuple of ints")
raise TypeError("Could not convert axis-input to tuple of ints")
if length is not None:
# shift negative indices to positive ones
......
......@@ -49,13 +49,14 @@ class RGRGTransformation(Transformation):
The axes along which the transformation should take place
"""
fct=1.
if self._transform.codomain.harmonic:
# correct for forward fft.
# naively one would set power to 0.5 here in order to
# apply effectively a factor of 1/sqrt(N) to the field.
# BUT: the pixel volumes of the domain and codomain are different.
# Hence, in order to produce the same scalar product, power===1.
val = self._transform.domain.weight(val, power=1, axes=axes)
fct *= self._transform.domain.weight()
# Perform the transformation
if issubclass(val.dtype.type, np.complexfloating):
......@@ -80,6 +81,7 @@ class RGRGTransformation(Transformation):
if not self._transform.codomain.harmonic:
# correct for inverse fft.
# See discussion above.
Tval = self._transform.codomain.weight(Tval, power=-1, axes=axes)
fct /= self._transform.codomain.weight()
Tval *= fct
return Tval
......@@ -6,6 +6,7 @@ import numpy as np
from ..endomorphic_operator import EndomorphicOperator
from ..fft_operator import FFTOperator
# MR FIXME: big optimization potential by precomputing kernel in constructor!
class FFTSmoothingOperator(EndomorphicOperator):
......
......@@ -109,31 +109,12 @@ class GLSpace(Space):
return 4 * np.pi
def copy(self):
return self.__class__(nlat=self.nlat,
nlon=self.nlon)
return self.__class__(nlat=self.nlat, nlon=self.nlon)
def weight(self, x, power=1, axes=None, inplace=False):
def weight(self):
from pyHealpix import GL_weights
nlon = self.nlon
nlat = self.nlat
vol = GL_weights(nlat, nlon) ** np.float(power)
weight = np.array(list(itertools.chain.from_iterable(
itertools.repeat(x, nlon) for x in vol)))
if axes is not None:
# reshape the weight array to match the input shape
new_shape = np.ones(len(x.shape), dtype=np.int)
# we know len(axes) is always 1
new_shape[axes[0]] = len(weight)
weight = weight.reshape(new_shape)
if inplace:
x *= weight
result_x = x
else:
result_x = x * weight
return result_x
vol = GL_weights(self.nlat, self.nlon)
return np.outer(vol, np.ones(self.nlon,dtype=np.float64)).flatten()
# ---Added properties and methods---
......@@ -142,21 +123,18 @@ class GLSpace(Space):
""" Number of latitudinal bins (or rings) that are used for this
pixelization.
"""
return self._nlat
@property
def nlon(self):
""" Number of longitudinal bins that are used for this pixelization.
"""
return self._nlon
def _parse_nlat(self, nlat):
nlat = int(nlat)
if nlat < 1:
raise ValueError(
"nlat must be a positive number.")
raise ValueError("nlat must be a positive number.")
return nlat
def _parse_nlon(self, nlon):
......
......@@ -80,7 +80,6 @@ class HPSpace(Space):
def __init__(self, nside):
super(HPSpace, self).__init__()
self._nside = self._parse_nside(nside)
# ---Mandatory properties and methods---
......@@ -94,7 +93,7 @@ class HPSpace(Space):
@property
def shape(self):
return (np.int(12 * self.nside * self.nside),)
return (self.dim,)
@property
def dim(self):
......@@ -107,17 +106,8 @@ class HPSpace(Space):
def copy(self):
return self.__class__(nside=self.nside)
def weight(self, x, power=1, axes=None, inplace=False):
weight = ((4*np.pi) / (12*self.nside*self.nside)) ** np.float(power)
if inplace:
x *= weight
result_x = x
else:
result_x = x * weight
return result_x
def weight(self):
return np.pi / (3*self._nside*self._nside)
# ---Added properties and methods---
......
......@@ -116,11 +116,8 @@ class LMSpace(Space):
def copy(self):
return self.__class__(lmax=self.lmax)
def weight(self, x, power=1, axes=None, inplace=False):
if inplace:
return x
else:
return x.copy()
def weight(self):
return 1.
def get_distance_array(self):
lmax = self.lmax
......
......@@ -178,22 +178,8 @@ class PowerSpace(Space):
return self.__class__(harmonic_partner=self.harmonic_partner,
binbounds=self._binbounds)
def weight(self, x, power, axes, inplace=False):
reshaper = [1, ] * len(x.shape)
# we know len(axes) is always 1
reshaper[axes[0]] = self.shape[0]
weight = self.rho.reshape(reshaper)
if power != 1:
weight = weight ** np.float(power)
if inplace:
x *= weight
result_x = x
else:
result_x = x*weight
return result_x
def weight(self):
return 1.
def get_distance_array(self):
return self.kindex.copy()
......
......@@ -115,14 +115,8 @@ class RGSpace(Space):
distances=self.distances,
harmonic=self.harmonic)
def weight(self, x, power=1, axes=None, inplace=False):
weight = reduce(lambda x, y: x*y, self.distances) ** np.float(power)
if inplace:
x *= weight
result_x = x
else:
result_x = x*weight
return result_x
def weight(self):
return reduce(lambda x, y: x*y, self.distances)
def get_distance_array(self):
""" Calculates an n-dimensional array with its entries being the
......
......@@ -47,18 +47,8 @@ def get_weight_configs():
# for GLSpace(nlat=2, nlon=3)
weight_0 = np.array(list(itertools.chain.from_iterable(
itertools.repeat(x, 3) for x in wgt)))
w_0_x = np.random.rand(6)
w_0_res = w_0_x * weight_0
weight_1 = np.array(list(itertools.chain.from_iterable(
itertools.repeat(x, 3) for x in wgt)))
weight_1 = weight_1.reshape([1, 1, 6])
w_1_x = np.random.rand(32, 16, 6)
w_1_res = w_1_x * weight_1
return [
[w_0_x, 1, None, False, w_0_res],
[w_0_x.copy(), 1, None, True, w_0_res],
[w_1_x.copy(), 1, (2,), True, w_1_res],
[1, weight_0],
]
......@@ -84,9 +74,5 @@ class GLSpaceFunctionalityTests(unittest.TestCase):
assert_equal(getattr(g, key), value)
@expand(get_weight_configs())
def test_weight(self, x, power, axes, inplace, expected):
g = GLSpace(2)
res = g.weight(x, power, axes, inplace)
assert_almost_equal(res, expected)
if inplace:
assert_(x is res)
def test_weight(self, power, expected):
assert_almost_equal(GLSpace(2).weight(), expected)
......@@ -54,20 +54,6 @@ CONSTRUCTOR_CONFIGS = [
]
def get_weight_configs():
np.random.seed(42)
# for HPSpace(nside=2)
w_0_x = np.random.rand(48)
w_0_res = w_0_x * ((4 * np.pi) / 48)
w_1_res = w_0_x * (((4 * np.pi) / 48)**2)
return [
[w_0_x, 1, None, False, w_0_res],
[w_0_x.copy(), 1, None, True, w_0_res],
[w_0_x, 2, None, False, w_1_res],
]
class HPSpaceInterfaceTests(unittest.TestCase):
@expand([['nside', int]])
def test_property_ret_type(self, attribute, expected_type):
......@@ -86,10 +72,5 @@ class HPSpaceFunctionalityTests(unittest.TestCase):
for key, value in expected.items():
assert_equal(getattr(h, key), value)
@expand(get_weight_configs())
def test_weight(self, x, power, axes, inplace, expected):
h = HPSpace(2)
res = h.weight(x, power, axes, inplace)
assert_almost_equal(res, expected)
if inplace:
assert_(x is res)
def test_weight(self):
assert_almost_equal(HPSpace(2).weight(), np.pi/12)
......@@ -69,15 +69,6 @@ def get_distance_array_configs():
return [[5, da_0]]
def get_weight_configs():
np.random.seed(42)
w_0_x = np.random.rand(32, 16, 6)
return [
[w_0_x, 1, None, False, w_0_x],
[w_0_x.copy(), 1, None, True, w_0_x]
]
class LMSpaceInterfaceTests(unittest.TestCase):
@expand([['lmax', int],
['mmax', int],
......@@ -98,13 +89,8 @@ class LMSpaceFunctionalityTests(unittest.TestCase):
for key, value in expected.items():
assert_equal(getattr(l, key), value)
@expand(get_weight_configs())
def test_weight(self, x, power, axes, inplace, expected):
l = LMSpace(5)
res = l.weight(x, power, axes, inplace)
assert_almost_equal(res, expected)
if inplace:
assert_(x is res)
def test_weight(self):
assert_almost_equal(LMSpace(5).weight(), 1.)
@expand(get_distance_array_configs())
def test_distance_array(self, lmax, expected):
......
......@@ -80,24 +80,6 @@ def get_distance_array_configs():
]
def get_weight_configs():
np.random.seed(42)
# power 1
w_0_x = np.random.rand(32, 16, 6)
# RGSpace((4, 4), harmonic=True)
# using rho directly
weight_0 = np.array([1, 4, 4, 2, 4, 1])
weight_0 = weight_0.reshape([1, 1, 6])
w_0_res = w_0_x * weight_0
return [
[RGSpace((4, 4), harmonic=True),
w_0_x, 1, (2,), False, w_0_res],
[RGSpace((4, 4), harmonic=True),
w_0_x.copy(), 1, (2,), True, w_0_res],
]
class PowerSpaceInterfaceTest(unittest.TestCase):
@expand([
['harmonic_partner', Space],
......@@ -148,11 +130,6 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
p = PowerSpace(harmonic_partner=harmonic_partner)
assert_almost_equal(p.get_distance_array(), expected)
@expand(get_weight_configs())
def test_weight(self, harmonic_partner, x, power, axes,
inplace, expected):
p = PowerSpace(harmonic_partner=harmonic_partner)
res = p.weight(x, power, axes, inplace)
assert_almost_equal(res, expected)
if inplace:
assert_(x is res)
def test_weight(self):
p = PowerSpace(harmonic_partner=RGSpace(10,harmonic=True))
assert_almost_equal(p.weight(),1.)
......@@ -90,25 +90,15 @@ def get_distance_array_configs():
def get_weight_configs():
np.random.seed(42)
# power 1
w_0_x = np.random.rand(32, 12, 6)
# for RGSpace(shape=(11,11), distances=None, harmonic=False)
w_0_res = w_0_x * (1/11 * 1/11)
# for RGSpace(shape=(11, 11), distances=(1.3,1.3), harmonic=False)
w_1_res = w_0_x * (1.3 * 1.3)
# for RGSpace(shape=(11,11), distances=None, harmonic=True)
w_2_res = w_0_x * (1.0 * 1.0)
# for RGSpace(shape=(11,11), distances=(1.3, 1,3), harmonic=True)
w_3_res = w_0_x * (1.3 * 1.3)
return [
[(11, 11), None, False, w_0_x, 1, None, False, w_0_res],
[(11, 11), None, False, w_0_x.copy(), 1, None, True, w_0_res],
[(11, 11), (1.3, 1.3), False, w_0_x, 1, None, False, w_1_res],
[(11, 11), (1.3, 1.3), False, w_0_x.copy(), 1, None, True, w_1_res],
[(11, 11), None, True, w_0_x, 1, None, False, w_2_res],
[(11, 11), None, True, w_0_x.copy(), 1, None, True, w_2_res],
[(11, 11), (1.3, 1.3), True, w_0_x, 1, None, False, w_3_res],
[(11, 11), (1.3, 1.3), True, w_0_x.copy(), 1, None, True, w_3_res]
[(11, 11), None, False, 1],
[(11, 11), None, False, 1],
[(11, 11), (1.3, 1.3), False, 1],
[(11, 11), (1.3, 1.3), False, 1],
[(11, 11), None, True, 1],
[(11, 11), None, True, 1],
[(11, 11), (1.3, 1.3), True, 1],
[(11, 11), (1.3, 1.3), True, 1]
]
......@@ -133,10 +123,6 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
assert_almost_equal(r.get_distance_array(), expected)
@expand(get_weight_configs())
def test_weight(self, shape, distances, harmonic, x, power, axes,
inplace, expected):
def test_weight(self, shape, distances, harmonic, power):
r = RGSpace(shape=shape, distances=distances, harmonic=harmonic)
res = r.weight(x, power, axes, inplace)
assert_almost_equal(res, expected)
if inplace:
assert_(x is res)
assert_almost_equal(r.weight(), np.prod(r.distances)**power)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment