Commit 643d308b authored by Martin Reinecke's avatar Martin Reinecke
Browse files

more dobj; kindex->k_lengths

parent 0674df96
...@@ -5,7 +5,7 @@ np.random.seed(42) ...@@ -5,7 +5,7 @@ np.random.seed(42)
#np.seterr(all="raise",under="ignore") #np.seterr(all="raise",under="ignore")
def plot_parameters(m, t, p, p_d): def plot_parameters(m, t, p, p_d):
x = np.log(t.domain[0].kindex) x = np.log(t.domain[0].k_lengths)
m = fft.adjoint_times(m) m = fft.adjoint_times(m)
t = t.val.real t = t.val.real
p = p.val.real p = p.val.real
...@@ -59,7 +59,7 @@ if __name__ == "__main__": ...@@ -59,7 +59,7 @@ if __name__ == "__main__":
S = ift.create_power_operator(h_space, power_spectrum=p_spec) S = ift.create_power_operator(h_space, power_spectrum=p_spec)
# Draw a sample sh from the prior distribution in harmonic space # Draw a sample sh from the prior distribution in harmonic space
sp = ift.Field(p_space, val=p_spec(p_space.kindex)) sp = ift.Field(p_space, val=p_spec(p_space.k_lengths))
sh = sp.power_synthesize(real_signal=True) sh = sp.power_synthesize(real_signal=True)
# Choose the measurement instrument # Choose the measurement instrument
...@@ -101,7 +101,7 @@ if __name__ == "__main__": ...@@ -101,7 +101,7 @@ if __name__ == "__main__":
def ps0(k): def ps0(k):
return (1./(1.+k)**2) return (1./(1.+k)**2)
t0 = ift.Field(p_space, val=np.log(1./(1+p_space.kindex)**2)) t0 = ift.Field(p_space, val=np.log(1./(1+p_space.k)_lengths)**2))
for i in range(500): for i in range(500):
S0 = ift.create_power_operator(h_space, power_spectrum=ps0) S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
......
...@@ -26,7 +26,7 @@ if __name__ == "__main__": ...@@ -26,7 +26,7 @@ if __name__ == "__main__":
# Creating the mock signal |\label{code:wf_mock_signal}| # Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum) S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.Field(power_space, val=power_spectrum(power_space.kindex)) mock_power = ift.Field(power_space, val=power_spectrum(power_space.k_lengths))
mock_signal = fft(mock_power.power_synthesize(real_signal=True)) mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response # Setting up an exemplary response
......
...@@ -25,7 +25,7 @@ if __name__ == "__main__": ...@@ -25,7 +25,7 @@ if __name__ == "__main__":
fft_1 = ift.FFTOperator(harmonic_space_1, target=signal_space_1) fft_1 = ift.FFTOperator(harmonic_space_1, target=signal_space_1)
power_space_1 = ift.PowerSpace(harmonic_space_1) power_space_1 = ift.PowerSpace(harmonic_space_1)
mock_power_1 = ift.Field(power_space_1, val=power_spectrum_1(power_space_1.kindex)) mock_power_1 = ift.Field(power_space_1, val=power_spectrum_1(power_space_1.k_lengths))
...@@ -48,7 +48,7 @@ if __name__ == "__main__": ...@@ -48,7 +48,7 @@ if __name__ == "__main__":
fft_2 = ift.FFTOperator(harmonic_space_2, target=signal_space_2) fft_2 = ift.FFTOperator(harmonic_space_2, target=signal_space_2)
power_space_2 = ift.PowerSpace(harmonic_space_2) power_space_2 = ift.PowerSpace(harmonic_space_2)
mock_power_2 = ift.Field(power_space_2, val=power_spectrum_2(power_space_2.kindex)) mock_power_2 = ift.Field(power_space_2, val=power_spectrum_2(power_space_2.k_lengths))
fft = ift.ComposedOperator((fft_1, fft_2)) fft = ift.ComposedOperator((fft_1, fft_2))
......
...@@ -25,7 +25,7 @@ if __name__ == "__main__": ...@@ -25,7 +25,7 @@ if __name__ == "__main__":
# Creating the mock signal |\label{code:wf_mock_signal}| # Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum) S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.Field(power_space, val=power_spectrum(power_space.kindex)) mock_power = ift.Field(power_space, val=power_spectrum(power_space.k_lengths))
mock_signal = fft(mock_power.power_synthesize(real_signal=True)) mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response # Setting up an exemplary response
......
...@@ -38,7 +38,7 @@ if __name__ == "__main__": ...@@ -38,7 +38,7 @@ if __name__ == "__main__":
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum) S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
np.random.seed(43) np.random.seed(43)
mock_power = ift.Field(power_space, val=power_spectrum(power_space.kindex)) mock_power = ift.Field(power_space, val=power_spectrum(power_space.k_lengths))
mock_harmonic = mock_power.power_synthesize(real_signal=True) mock_harmonic = mock_power.power_synthesize(real_signal=True)
mock_harmonic = mock_harmonic.real mock_harmonic = mock_harmonic.real
mock_signal = fft(mock_harmonic) mock_signal = fft(mock_harmonic)
......
...@@ -50,7 +50,7 @@ if __name__ == "__main__": ...@@ -50,7 +50,7 @@ if __name__ == "__main__":
S = ift.create_power_operator(h_space, power_spectrum=p_spec) S = ift.create_power_operator(h_space, power_spectrum=p_spec)
# Drawing a sample sh from the prior distribution in harmonic space # Drawing a sample sh from the prior distribution in harmonic space
sp = ift.Field(p_space, val=p_spec(p_space.kindex)) sp = ift.Field(p_space, val=p_spec(p_space.k_lengths))
sh = sp.power_synthesize(real_signal=True) sh = sp.power_synthesize(real_signal=True)
ss = fft.adjoint_times(sh) ss = fft.adjoint_times(sh)
......
...@@ -19,6 +19,7 @@ ...@@ -19,6 +19,7 @@
from __future__ import division from __future__ import division
import numpy as np import numpy as np
from .field import Field from .field import Field
from . import dobj
__all__ = ['cos', 'sin', 'cosh', 'sinh', 'tan', 'tanh', 'arccos', 'arcsin', __all__ = ['cos', 'sin', 'cosh', 'sinh', 'tan', 'tanh', 'arccos', 'arcsin',
...@@ -39,68 +40,68 @@ def _math_helper(x, function, out): ...@@ -39,68 +40,68 @@ def _math_helper(x, function, out):
def cos(x, out=None): def cos(x, out=None):
return _math_helper(x, np.cos, out) return _math_helper(x, dobj.cos, out)
def sin(x, out=None): def sin(x, out=None):
return _math_helper(x, np.sin, out) return _math_helper(x, dobj.sin, out)
def cosh(x, out=None): def cosh(x, out=None):
return _math_helper(x, np.cosh, out) return _math_helper(x, dobj.cosh, out)
def sinh(x, out=None): def sinh(x, out=None):
return _math_helper(x, np.sinh, out) return _math_helper(x, dobj.sinh, out)
def tan(x, out=None): def tan(x, out=None):
return _math_helper(x, np.tan, out) return _math_helper(x, dobj.tan, out)
def tanh(x, out=None): def tanh(x, out=None):
return _math_helper(x, np.tanh, out) return _math_helper(x, dobj.tanh, out)
def arccos(x, out=None): def arccos(x, out=None):
return _math_helper(x, np.arccos, out) return _math_helper(x, dobj.arccos, out)
def arcsin(x, out=None): def arcsin(x, out=None):
return _math_helper(x, np.arcsin, out) return _math_helper(x, dobj.arcsin, out)
def arccosh(x, out=None): def arccosh(x, out=None):
return _math_helper(x, np.arccosh, out) return _math_helper(x, dobj.arccosh, out)
def arcsinh(x, out=None): def arcsinh(x, out=None):
return _math_helper(x, np.arcsinh, out) return _math_helper(x, dobj.arcsinh, out)
def arctan(x, out=None): def arctan(x, out=None):
return _math_helper(x, np.arctan, out) return _math_helper(x, dobj.arctan, out)
def arctanh(x, out=None): def arctanh(x, out=None):
return _math_helper(x, np.arctanh, out) return _math_helper(x, dobj.arctanh, out)
def sqrt(x, out=None): def sqrt(x, out=None):
return _math_helper(x, np.sqrt, out) return _math_helper(x, dobj.sqrt, out)
def exp(x, out=None): def exp(x, out=None):
return _math_helper(x, np.exp, out) return _math_helper(x, dobj.exp, out)
def log(x, out=None): def log(x, out=None):
return _math_helper(x, np.log, out) return _math_helper(x, dobj.log, out)
def conjugate(x, out=None): def conjugate(x, out=None):
return _math_helper(x, np.conjugate, out) return _math_helper(x, dobj.conjugate, out)
def conj(x, out=None): def conj(x, out=None):
return _math_helper(x, np.conj, out) return _math_helper(x, dobj.conj, out)
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
import numpy as np import numpy as np
from numpy import ndarray as data_object from numpy import ndarray as data_object
from numpy import full, empty, sqrt, ones, zeros, vdot, abs from numpy import full, empty, sqrt, ones, zeros, vdot, abs, bincount
def from_object(object, dtype=None, copy=True): def from_object(object, dtype=None, copy=True):
return np.array(object, dtype=dtype, copy=copy) return np.array(object, dtype=dtype, copy=copy)
...@@ -50,7 +50,7 @@ class LaplaceOperator(EndomorphicOperator): ...@@ -50,7 +50,7 @@ class LaplaceOperator(EndomorphicOperator):
self._logarithmic = bool(logarithmic) self._logarithmic = bool(logarithmic)
pos = self.domain[0].kindex.copy() pos = self.domain[0].k_lengths.copy()
if self.logarithmic: if self.logarithmic:
pos[1:] = np.log(pos[1:]) pos[1:] = np.log(pos[1:])
pos[0] = pos[1]-1. pos[0] = pos[1]-1.
......
...@@ -197,7 +197,7 @@ def plot(f, **kwargs): ...@@ -197,7 +197,7 @@ def plot(f, **kwargs):
_makeplot(kwargs.get("name")) _makeplot(kwargs.get("name"))
return return
elif isinstance(dom, PowerSpace): elif isinstance(dom, PowerSpace):
xcoord = dom.kindex xcoord = dom.k_lengths
ycoord = f.val ycoord = f.val
plt.xscale('log') plt.xscale('log')
plt.yscale('log') plt.yscale('log')
......
...@@ -20,6 +20,7 @@ import numpy as np ...@@ -20,6 +20,7 @@ import numpy as np
from ...spaces.space import Space from ...spaces.space import Space
from functools import reduce from functools import reduce
from ... import dobj
class PowerSpace(Space): class PowerSpace(Space):
...@@ -45,10 +46,10 @@ class PowerSpace(Space): ...@@ -45,10 +46,10 @@ class PowerSpace(Space):
Attributes Attributes
---------- ----------
pindex : numpy.ndarray pindex : data object
This holds the information which pixel of the harmonic partner gets This holds the information which pixel of the harmonic partner gets
mapped to which power bin mapped to which power bin
kindex : numpy.ndarray k_lengths : numpy.ndarray
Sorted array of all k-modes. Sorted array of all k-modes.
rho : numpy.ndarray rho : numpy.ndarray
The amount of k-modes that get mapped to one power bin is given by The amount of k-modes that get mapped to one power bin is given by
...@@ -161,17 +162,17 @@ class PowerSpace(Space): ...@@ -161,17 +162,17 @@ 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 np.any(temp_rho == 0), "empty bins detected"
temp_kindex = np.bincount(temp_pindex.ravel(), temp_k_lengths = np.bincount(temp_pindex.ravel(),
weights=k_length_array.ravel()) \ weights=k_length_array.ravel()) \
/ temp_rho / temp_rho
self._powerIndexCache[key] = (binbounds, self._powerIndexCache[key] = (binbounds,
temp_pindex, temp_pindex,
temp_kindex, temp_k_lengths,
temp_rho) temp_rho)
(self._binbounds, self._pindex, self._kindex, self._rho) = \ (self._binbounds, self._pindex, self._k_lengths, self._rho) = \
self._powerIndexCache[key] self._powerIndexCache[key]
@staticmethod @staticmethod
...@@ -193,7 +194,7 @@ class PowerSpace(Space): ...@@ -193,7 +194,7 @@ class PowerSpace(Space):
@property @property
def shape(self): def shape(self):
return self.kindex.shape return self.k_lengths.shape
@property @property
def dim(self): def dim(self):
...@@ -216,7 +217,7 @@ class PowerSpace(Space): ...@@ -216,7 +217,7 @@ class PowerSpace(Space):
return np.asarray(self.rho, dtype=np.float64) return np.asarray(self.rho, dtype=np.float64)
def get_k_length_array(self): def get_k_length_array(self):
return self.kindex.copy() return self.k_lengths.copy()
def get_fft_smoothing_kernel_function(self, sigma): def get_fft_smoothing_kernel_function(self, sigma):
raise NotImplementedError( raise NotImplementedError(
...@@ -236,16 +237,16 @@ class PowerSpace(Space): ...@@ -236,16 +237,16 @@ class PowerSpace(Space):
@property @property
def pindex(self): def pindex(self):
""" A numpy.ndarray having the shape of the harmonic partner """ A data object having the shape of the harmonic partner
space containing the indices of the power bin a pixel belongs to. space containing the indices of the power bin a pixel belongs to.
""" """
return self._pindex return self._pindex
@property @property
def kindex(self): def k_lengths(self):
""" Sorted array of all k-modes. """ Sorted array of all k-modes.
""" """
return self._kindex return self._k_lengths
@property @property
def rho(self): def rho(self):
......
...@@ -58,7 +58,7 @@ def create_power_operator(domain, power_spectrum, dtype=None): ...@@ -58,7 +58,7 @@ def create_power_operator(domain, power_spectrum, dtype=None):
raise TypeError("power_spectrum must be callable") raise TypeError("power_spectrum must be callable")
power_domain = PowerSpace(domain) power_domain = PowerSpace(domain)
fp = Field(power_domain, val=power_spectrum(power_domain.kindex), fp = Field(power_domain, val=power_spectrum(power_domain.k_lengths),
dtype=dtype) dtype=dtype)
f = fp.power_synthesize_special() f = fp.power_synthesize_special()
......
...@@ -63,11 +63,11 @@ class Test_Functionality(unittest.TestCase): ...@@ -63,11 +63,11 @@ class Test_Functionality(unittest.TestCase):
p1 = PowerSpace(space1) p1 = PowerSpace(space1)
spec1 = lambda k: 42/(1+k)**2 spec1 = lambda k: 42/(1+k)**2
fp1 = Field(p1, val=spec1(p1.kindex)) fp1 = Field(p1, val=spec1(p1.k_lengths))
p2 = PowerSpace(space2) p2 = PowerSpace(space2)
spec2 = lambda k: 42/(1+k)**3 spec2 = lambda k: 42/(1+k)**3
fp2 = Field(p2, val=spec2(p2.kindex)) fp2 = Field(p2, val=spec2(p2.k_lengths))
outer = np.outer(fp1.val, fp2.val) outer = np.outer(fp1.val, fp2.val)
fp = Field((p1, p2), val=outer) fp = Field((p1, p2), val=outer)
......
...@@ -56,7 +56,7 @@ CONSTRUCTOR_CONFIGS = [ ...@@ -56,7 +56,7 @@ CONSTRUCTOR_CONFIGS = [
'harmonic_partner': RGSpace((8,), harmonic=True), 'harmonic_partner': RGSpace((8,), harmonic=True),
'binbounds': None, 'binbounds': None,
'pindex': np.array([0, 1, 2, 3, 4, 3, 2, 1]), 'pindex': np.array([0, 1, 2, 3, 4, 3, 2, 1]),
'kindex': np.array([0., 1., 2., 3., 4.]), 'k_lengths': np.array([0., 1., 2., 3., 4.]),
'rho': np.array([1, 2, 2, 2, 1]), 'rho': np.array([1, 2, 2, 2, 1]),
}], }],
[RGSpace((8,), harmonic=True), True, None, None, { [RGSpace((8,), harmonic=True), True, None, None, {
...@@ -67,7 +67,7 @@ CONSTRUCTOR_CONFIGS = [ ...@@ -67,7 +67,7 @@ CONSTRUCTOR_CONFIGS = [
'harmonic_partner': RGSpace((8,), harmonic=True), 'harmonic_partner': RGSpace((8,), harmonic=True),
'binbounds': (0.5, 1.3228756555322954, 3.5), 'binbounds': (0.5, 1.3228756555322954, 3.5),
'pindex': np.array([0, 1, 2, 2, 3, 2, 2, 1]), 'pindex': np.array([0, 1, 2, 2, 3, 2, 2, 1]),
'kindex': np.array([0., 1., 2.5, 4.]), 'k_lengths': np.array([0., 1., 2.5, 4.]),
'rho': np.array([1, 2, 4, 1]), 'rho': np.array([1, 2, 4, 1]),
}], }],
] ]
...@@ -85,7 +85,7 @@ class PowerSpaceInterfaceTest(unittest.TestCase): ...@@ -85,7 +85,7 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
['harmonic_partner', Space], ['harmonic_partner', Space],
['binbounds', type(None)], ['binbounds', type(None)],
['pindex', np.ndarray], ['pindex', np.ndarray],
['kindex', np.ndarray], ['k_lengths', np.ndarray],
['rho', np.ndarray], ['rho', np.ndarray],
]) ])
def test_property_ret_type(self, attribute, expected_type): def test_property_ret_type(self, attribute, expected_type):
......
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