Commit a13b4c77 authored by Martin Reinecke's avatar Martin Reinecke

work in progress

parent 550a546a
Pipeline #20660 failed with stage
in 3 minutes and 54 seconds
import numpy as np
from .random import Random
class data_object(object):
def __init__(self, npdata):
self._data = np.asarray(npdata)
def __getitem__(self, key):
res = self._data[key]
return res if np.isscalar(res) else data_object(res)
def __setitem__(self, key, value):
self._data[key] = value
@property
def dtype(self):
return self._data.dtype
@property
def shape(self):
return self._data.shape
@property
def size(self):
return self._data.size
@property
def real(self):
return data_object(self._data.real)
@property
def imag(self):
return data_object(self._data.imag)
def _contraction_helper(self, op, axis):
if axis is None:
return getattr(self._data, op)()
# perform the contraction on the data
data = getattr(self._data, op)(axis=axis)
# check if the result is scalar or if a result_field must be constr.
if np.isscalar(data):
return data
else:
return data_object(data)
def sum(self, axis=None):
return self._contraction_helper("sum", axis)
def _binary_helper(self, other, op):
a=self._data
if isinstance(other, data_object):
b=other._data
#if a.shape != b.shape:
# print("shapes are incompatible.")
else:
b=other
tval = getattr(a, op)(b)
return self if tval is a else data_object(tval)
def __add__(self, other):
return self._binary_helper(other, op='__add__')
def __radd__(self, other):
return self._binary_helper(other, op='__radd__')
def __iadd__(self, other):
return self._binary_helper(other, op='__iadd__')
def __sub__(self, other):
return self._binary_helper(other, op='__sub__')
def __rsub__(self, other):
return self._binary_helper(other, op='__rsub__')
def __isub__(self, other):
return self._binary_helper(other, op='__isub__')
def __mul__(self, other):
return self._binary_helper(other, op='__mul__')
def __rmul__(self, other):
return self._binary_helper(other, op='__rmul__')
def __imul__(self, other):
return self._binary_helper(other, op='__imul__')
def __div__(self, other):
return self._binary_helper(other, op='__div__')
def __rdiv__(self, other):
return self._binary_helper(other, op='__rdiv__')
def __truediv__(self, other):
return self._binary_helper(other, op='__truediv__')
def __rtruediv__(self, other):
return self._binary_helper(other, op='__rtruediv__')
def __pow__(self, other):
return self._binary_helper(other, op='__pow__')
def __rpow__(self, other):
return self._binary_helper(other, op='__rpow__')
def __ipow__(self, other):
return self._binary_helper(other, op='__ipow__')
def __eq__(self, other):
return self._binary_helper(other, op='__eq__')
def __ne__(self, other):
return self._binary_helper(other, op='__ne__')
def __neg__(self):
return data_object(-self._data)
def __abs__(self):
return data_object(np.abs(self._data))
def ravel(self):
return data_object(self._data.ravel())
def reshape(self, shape):
return data_object(self._data.reshape(shape))
def all(self):
return self._data.all()
def any(self):
return self._data.any()
def full(shape, fill_value, dtype=None):
return data_object(np.full(shape, fill_value, dtype))
def empty(shape, dtype=np.float):
return data_object(np.empty(shape, dtype))
def zeros(shape, dtype=np.float):
return data_object(np.zeros(shape, dtype))
def ones(shape, dtype=np.float):
return data_object(np.ones(shape, dtype))
def empty_like(a, dtype=None):
return data_object(np.empty_like(a._data, dtype))
def vdot(a,b):
return np.vdot(a._data, b._data)
def abs(a, out=None):
if out is None:
out = empty_like(a)
np.abs(a._data, out=out._data)
return out
def exp(a, out=None):
if out is None:
out = empty_like(a)
np.exp(a._data, out=out._data)
return out
def log(a, out=None):
if out is None:
out = empty_like(a)
np.log(a._data, out=out._data)
return out
def sqrt(a, out=None):
if out is None:
out = empty_like(a)
np.sqrt(a._data, out=out._data)
return out
def bincount(x, weights=None, minlength=0):
if weights is not None:
weights = weights._data
res = np.bincount(x._data, weights, minlength)
return data_object(res)
def from_object(object, dtype=None, copy=True):
return data_object(np.array(object._data, dtype=dtype, copy=copy))
def from_random(random_type, shape, dtype=np.float64, **kwargs):
generator_function = getattr(Random, random_type)
return data_object(generator_function(dtype=dtype, shape=shape, **kwargs))
def to_ndarray(arr):
return arr._data
def from_ndarray(arr):
return data_object(arr)
...@@ -13,3 +13,9 @@ def from_object(object, dtype=None, copy=True): ...@@ -13,3 +13,9 @@ def from_object(object, dtype=None, copy=True):
def from_random(random_type, shape, dtype=np.float64, **kwargs): def from_random(random_type, shape, dtype=np.float64, **kwargs):
generator_function = getattr(Random, random_type) generator_function = getattr(Random, random_type)
return generator_function(dtype=dtype, shape=shape, **kwargs) return generator_function(dtype=dtype, shape=shape, **kwargs)
def to_ndarray(arr):
return arr
def from_ndarray(arr):
return np.asarray(arr)
from .data_objects.numpy_do import * from .data_objects.my_own_do import *
...@@ -308,7 +308,7 @@ class Field(object): ...@@ -308,7 +308,7 @@ class Field(object):
if np.isscalar(wgt): if np.isscalar(wgt):
fct *= wgt fct *= wgt
else: else:
new_shape = dobj.ones(len(self.shape), dtype=np.int) new_shape = np.ones(len(self.shape), dtype=np.int)
new_shape[self.domain.axes[ind][0]: new_shape[self.domain.axes[ind][0]:
self.domain.axes[ind][-1]+1] = wgt.shape self.domain.axes[ind][-1]+1] = wgt.shape
wgt = wgt.reshape(new_shape) wgt = wgt.reshape(new_shape)
......
...@@ -23,6 +23,7 @@ from ..field import Field ...@@ -23,6 +23,7 @@ from ..field import Field
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from .endomorphic_operator import EndomorphicOperator from .endomorphic_operator import EndomorphicOperator
from ..nifty_utilities import cast_iseq_to_tuple from ..nifty_utilities import cast_iseq_to_tuple
from .. import dobj
class DiagonalOperator(EndomorphicOperator): class DiagonalOperator(EndomorphicOperator):
""" NIFTY class for diagonal operators. """ NIFTY class for diagonal operators.
...@@ -146,7 +147,7 @@ class DiagonalOperator(EndomorphicOperator): ...@@ -146,7 +147,7 @@ class DiagonalOperator(EndomorphicOperator):
reshaper = [shp if i in active_axes else 1 reshaper = [shp if i in active_axes else 1
for i, shp in enumerate(x.shape)] for i, shp in enumerate(x.shape)]
reshaped_local_diagonal = np.reshape(self._diagonal.val, reshaper) reshaped_local_diagonal = self._diagonal.val.reshape(reshaper)
# here the actual multiplication takes place # here the actual multiplication takes place
return Field(x.domain, val=operation(reshaped_local_diagonal)(x.val)) return Field(x.domain, val=operation(reshaped_local_diagonal)(x.val))
...@@ -20,6 +20,7 @@ from __future__ import division ...@@ -20,6 +20,7 @@ from __future__ import division
import numpy as np import numpy as np
from .. import nifty_utilities as utilities from .. import nifty_utilities as utilities
from ..low_level_library import hartley from ..low_level_library import hartley
from ..dobj import to_ndarray as to_np, from_ndarray as from_np
class Transformation(object): class Transformation(object):
def __init__(self, domain, codomain): def __init__(self, domain, codomain):
...@@ -68,10 +69,10 @@ class RGRGTransformation(Transformation): ...@@ -68,10 +69,10 @@ class RGRGTransformation(Transformation):
# Perform the transformation # Perform the transformation
if issubclass(val.dtype.type, np.complexfloating): if issubclass(val.dtype.type, np.complexfloating):
Tval = hartley(val.real, axes) \ Tval = from_np(hartley(to_np(val.real), axes) \
+ 1j*hartley(val.imag, axes) + 1j*hartley(to_np(val.imag), axes))
else: else:
Tval = hartley(val, axes) Tval = from_np(hartley(to_np(val), axes))
return Tval, fct return Tval, fct
......
...@@ -57,7 +57,9 @@ class PowerProjectionOperator(LinearOperator): ...@@ -57,7 +57,9 @@ class PowerProjectionOperator(LinearOperator):
x.domain.collapsed_shape_for_domain(self._space)) x.domain.collapsed_shape_for_domain(self._space))
out = dobj.zeros(self._target.collapsed_shape_for_domain(self._space), out = dobj.zeros(self._target.collapsed_shape_for_domain(self._space),
dtype=x.dtype) dtype=x.dtype)
np.add.at(out, (slice(None), pindex.ravel(), slice(None)), arr) out = dobj.to_ndarray(out)
np.add.at(out, (slice(None), dobj.to_ndarray(pindex.ravel()), slice(None)), dobj.to_ndarray(arr))
out = dobj.from_ndarray(out)
return Field(self._target, out.reshape(self._target.shape))\ return Field(self._target, out.reshape(self._target.shape))\
.weight(-1, spaces=self._space) .weight(-1, spaces=self._space)
...@@ -65,7 +67,7 @@ class PowerProjectionOperator(LinearOperator): ...@@ -65,7 +67,7 @@ class PowerProjectionOperator(LinearOperator):
pindex = self._target[self._space].pindex pindex = self._target[self._space].pindex
pindex = pindex.reshape((1, pindex.size, 1)) pindex = pindex.reshape((1, pindex.size, 1))
arr = x.val.reshape(x.domain.collapsed_shape_for_domain(self._space)) arr = x.val.reshape(x.domain.collapsed_shape_for_domain(self._space))
out = arr[(slice(None), pindex.ravel(), slice(None))] out = arr[(slice(None), pindex.ravel()._data, slice(None))]
return Field(self._domain, out.reshape(self._domain.shape)) return Field(self._domain, out.reshape(self._domain.shape))
@property @property
......
...@@ -141,10 +141,10 @@ class PowerSpace(Space): ...@@ -141,10 +141,10 @@ class PowerSpace(Space):
harmonic_partner=self.harmonic_partner, harmonic_partner=self.harmonic_partner,
k_length_array=k_length_array, k_length_array=k_length_array,
binbounds=binbounds) binbounds=binbounds)
temp_rho = np.bincount(temp_pindex.ravel()) temp_rho = dobj.bincount(temp_pindex.ravel())
assert not np.any(temp_rho == 0), "empty bins detected" assert not (temp_rho == 0).any(), "empty bins detected"
temp_k_lengths = np.bincount(temp_pindex.ravel(), temp_k_lengths = dobj.bincount(temp_pindex.ravel(),
weights=k_length_array.ravel()) \ weights=dobj.from_ndarray(k_length_array.ravel())) \
/ temp_rho / temp_rho
temp_dvol = temp_rho*pdvol temp_dvol = temp_rho*pdvol
self._powerIndexCache[key] = (binbounds, self._powerIndexCache[key] = (binbounds,
...@@ -160,7 +160,7 @@ class PowerSpace(Space): ...@@ -160,7 +160,7 @@ class PowerSpace(Space):
if binbounds is None: if binbounds is None:
tmp = harmonic_partner.get_unique_k_lengths() tmp = harmonic_partner.get_unique_k_lengths()
binbounds = 0.5*(tmp[:-1]+tmp[1:]) binbounds = 0.5*(tmp[:-1]+tmp[1:])
return np.searchsorted(binbounds, k_length_array) return dobj.from_ndarray(np.searchsorted(binbounds, k_length_array))
# ---Mandatory properties and methods--- # ---Mandatory properties and methods---
......
...@@ -27,6 +27,7 @@ from . import Space,\ ...@@ -27,6 +27,7 @@ from . import Space,\
sqrt,\ sqrt,\
DomainTuple DomainTuple
from . import nifty_utilities as utilities from . import nifty_utilities as utilities
from . import dobj
__all__ = ['power_analyze', __all__ = ['power_analyze',
'power_synthesize', 'power_synthesize',
...@@ -222,7 +223,8 @@ def create_power_field(domain, power_spectrum, dtype=None): ...@@ -222,7 +223,8 @@ def create_power_field(domain, power_spectrum, dtype=None):
fp = Field(power_domain, val=power_spectrum.val, dtype=dtype) fp = Field(power_domain, val=power_spectrum.val, dtype=dtype)
else: else:
power_domain = PowerSpace(domain) power_domain = PowerSpace(domain)
fp = Field(power_domain, val=power_spectrum(power_domain.k_lengths), fp = Field(power_domain,
val=dobj.from_ndarray(power_spectrum(dobj.to_ndarray(power_domain.k_lengths))),
dtype=dtype) dtype=dtype)
P = PowerProjectionOperator(domain, power_domain) P = PowerProjectionOperator(domain, power_domain)
f = P.adjoint_times(fp) f = P.adjoint_times(fp)
...@@ -258,7 +260,7 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None): ...@@ -258,7 +260,7 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
""" """
domain = DomainTuple.make(domain) domain = DomainTuple.make(domain)
if space is None: if space is None:
if len(domain)!=1: if len(domain) != 1:
raise ValueError("space keyword must be set") raise ValueError("space keyword must be set")
else: else:
space = 0 space = 0
......
...@@ -17,111 +17,108 @@ ...@@ -17,111 +17,108 @@
# and financially supported by the Studienstiftung des deutschen Volkes. # and financially supported by the Studienstiftung des deutschen Volkes.
import unittest import unittest
import numpy as np import numpy as np
from numpy.testing import assert_equal,\ from numpy.testing import assert_equal, assert_allclose
assert_allclose
from itertools import product from itertools import product
from nifty2go import Field,\ import nifty2go as ift
RGSpace,\ from nifty2go.dobj import to_ndarray as to_np, from_ndarray as from_np
LMSpace,\
PowerSpace,\
DomainTuple,\
DiagonalOperator
from nifty2go.sugar import create_power_field, power_analyze, power_synthesize
from test.common import expand from test.common import expand
SPACES = [RGSpace((4,)), RGSpace((5))] SPACES = [ift.RGSpace((4,)), ift.RGSpace((5))]
SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES] SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES]
class Test_Interface(unittest.TestCase): class Test_Interface(unittest.TestCase):
@expand(product(SPACE_COMBINATIONS, @expand(product(SPACE_COMBINATIONS,
[['domain', DomainTuple], [['domain', ift.DomainTuple],
['val', np.ndarray], ['val', ift.dobj.data_object],
['shape', tuple], ['shape', tuple],
['dim', (np.int, np.int64)]])) ['dim', (np.int, np.int64)]]))
def test_return_types(self, domain, attribute_desired_type): def test_return_types(self, domain, attribute_desired_type):
attribute = attribute_desired_type[0] attribute = attribute_desired_type[0]
desired_type = attribute_desired_type[1] desired_type = attribute_desired_type[1]
f = Field(domain=domain) f = ift.Field(domain=domain)
assert_equal(isinstance(getattr(f, attribute), desired_type), True) assert_equal(isinstance(getattr(f, attribute), desired_type), True)
def _spec1(k):
return 42/(1.+k)**2
def _spec2(k):
return 42/(1.+k)**3
class Test_Functionality(unittest.TestCase): class Test_Functionality(unittest.TestCase):
@expand(product([RGSpace((8,), harmonic=True), @expand(product([ift.RGSpace((8,), harmonic=True),
RGSpace((8, 8), harmonic=True, distances=0.123)], ift.RGSpace((8, 8), harmonic=True, distances=0.123)],
[RGSpace((8,), harmonic=True), [ift.RGSpace((8,), harmonic=True),
LMSpace(12)])) ift.LMSpace(12)]))
def test_power_synthesize_analyze(self, space1, space2): def test_power_synthesize_analyze(self, space1, space2):
np.random.seed(11) np.random.seed(11)
p1 = PowerSpace(space1) p1 = ift.PowerSpace(space1)
spec1 = lambda k: 42/(1+k)**2 fp1 = ift.Field(p1, val=from_np(_spec1(to_np(p1.k_lengths))))
fp1 = Field(p1, val=spec1(p1.k_lengths))
p2 = PowerSpace(space2) p2 = ift.PowerSpace(space2)
spec2 = lambda k: 42/(1+k)**3 fp2 = ift.Field(p2, val=from_np(_spec2(to_np(p2.k_lengths))))
fp2 = Field(p2, val=spec2(p2.k_lengths))
outer = np.outer(fp1.val, fp2.val) outer = from_np(np.outer(to_np(fp1.val), to_np(fp2.val)))
fp = Field((p1, p2), val=outer) fp = ift.Field((p1, p2), val=outer)
samples = 500 samples = 500
ps1 = 0. ps1 = 0.
ps2 = 0. ps2 = 0.
for ii in range(samples): for ii in range(samples):
sk = power_synthesize(fp, spaces=(0, 1), real_signal=True) sk = ift.power_synthesize(fp, spaces=(0, 1), real_signal=True)
sp = power_analyze(sk, spaces=(0, 1), keep_phase_information=False) sp = ift.power_analyze(sk, spaces=(0, 1),
keep_phase_information=False)
ps1 += sp.sum(spaces=1)/fp2.sum() ps1 += sp.sum(spaces=1)/fp2.sum()
ps2 += sp.sum(spaces=0)/fp1.sum() ps2 += sp.sum(spaces=0)/fp1.sum()
assert_allclose(ps1.val/samples, fp1.val, rtol=0.2) assert_allclose(to_np(ps1.val/samples), to_np(fp1.val), rtol=0.2)
assert_allclose(ps2.val/samples, fp2.val, rtol=0.2) assert_allclose(to_np(ps2.val/samples), to_np(fp2.val), rtol=0.2)
@expand(product([RGSpace((8,), harmonic=True), @expand(product([ift.RGSpace((8,), harmonic=True),
RGSpace((8, 8), harmonic=True, distances=0.123)], ift.RGSpace((8, 8), harmonic=True, distances=0.123)],
[RGSpace((8,), harmonic=True), [ift.RGSpace((8,), harmonic=True),
LMSpace(12)])) ift.LMSpace(12)]))
def test_DiagonalOperator_power_analyze(self, space1, space2): def test_DiagonalOperator_power_analyze(self, space1, space2):
np.random.seed(11) np.random.seed(11)
fulldomain = DomainTuple.make((space1, space2)) fulldomain = ift.DomainTuple.make((space1, space2))
p1 = PowerSpace(space1) p1 = ift.PowerSpace(space1)
spec1 = lambda k: 42/(1+k)**2 fp1 = ift.Field(p1, val=from_np(_spec1(to_np(p1.k_lengths))))
fp1 = Field(p1, val=spec1(p1.k_lengths))
p2 = PowerSpace(space2) p2 = ift.PowerSpace(space2)
spec2 = lambda k: 42/(1+k)**3 fp2 = ift.Field(p2, val=from_np(_spec2(to_np(p2.k_lengths))))
fp2 = Field(p2, val=spec2(p2.k_lengths))
S_1 = create_power_field(space1, lambda x: np.sqrt(spec1(x))) S_1 = ift.create_power_field(space1, lambda x: np.sqrt(_spec1(x)))
S_1 = DiagonalOperator(S_1, domain=fulldomain, spaces=0) S_1 = ift.DiagonalOperator(S_1, domain=fulldomain, spaces=0)
S_2 = create_power_field(space2, lambda x: np.sqrt(spec2(x))) S_2 = ift.create_power_field(space2, lambda x: np.sqrt(_spec2(x)))
S_2 = DiagonalOperator(S_2, domain=fulldomain, spaces=1) S_2 = ift.DiagonalOperator(S_2, domain=fulldomain, spaces=1)
samples = 500 samples = 500
ps1 = 0. ps1 = 0.
ps2 = 0. ps2 = 0.
for ii in range(samples): for ii in range(samples):
rand_k = Field.from_random('normal', domain=fulldomain) rand_k = ift.Field.from_random('normal', domain=fulldomain)
sk = S_1.times(S_2.times(rand_k)) sk = S_1.times(S_2.times(rand_k))
sp = power_analyze(sk, spaces=(0, 1), keep_phase_information=False) sp = ift.power_analyze(sk, spaces=(0, 1),
keep_phase_information=False)
ps1 += sp.sum(spaces=1)/fp2.sum() ps1 += sp.sum(spaces=1)/fp2.sum()
ps2 += sp.sum(spaces=0)/fp1.sum() ps2 += sp.sum(spaces=0)/fp1.sum()
assert_allclose(ps1.val/samples, fp1.val, rtol=0.2) assert_allclose(to_np(ps1.val/samples), to_np(fp1.val), rtol=0.2)
assert_allclose(ps2.val/samples, fp2.val, rtol=0.2) assert_allclose(to_np(ps2.val/samples), to_np(fp2.val), rtol=0.2)
def test_vdot(self): def test_vdot(self):
s = RGSpace((10,)) s = ift.RGSpace((10,))
f1 = Field.from_random("normal", domain=s, dtype=np.complex128) f1 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
f2 = Field.from_random("normal", domain=s, dtype=np.complex128) f2 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
assert_allclose(f1.vdot(f2), f1.vdot(f2, spaces=0)) assert_allclose(f1.vdot(f2), f1.vdot(f2, spaces=0))
assert_allclose(f1.vdot(f2), np.conj(f2.vdot(f1))) assert_allclose(f1.vdot(f2), np.conj(f2.vdot(f1)))
import unittest import unittest
import numpy as np import numpy as np
from numpy.testing import assert_allclose from numpy.testing import assert_allclose
import nifty2go as ift import nifty2go as ift