Commit 6ddae70b authored by Philipp Frank's avatar Philipp Frank
Browse files

fix everything with slope operator

parent 54e22df4
......@@ -78,6 +78,7 @@ class LogRGSpace(StructuredDomain):
self._t_0, True)
def get_k_length_array(self):
# FIXME This is simply wrong! Definitely not the final version.
ib = dobj.ibegin_from_shape(self._shape)
res = np.arange(self.local_shape[0], dtype=np.float64) + ib[0]
res = np.minimum(res, self.shape[0]-res)*self.bindistances[0]
......
......@@ -50,14 +50,10 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
# Prepare q_array
q_array = np.zeros((dim,) + shape)
if dim == 1:
ks = domain.get_k_length_array().to_global_data()
q_array = np.array([ks])
else:
for i in range(dim):
ks = np.minimum(shape[i] - np.arange(shape[i]) +
1, np.arange(shape[i])) * dist[i]
q_array[i] += ks.reshape((1,)*i + (shape[i],) + (1,)*(dim-i-1))
for i in range(dim):
ks = np.zeros(shape[i])
ks[1:] = np.minimum(shape[i] - 1 - np.arange(shape[i]-1), np.arange(shape[i]-1)) * dist[i]
q_array[i] += ks.reshape((1,)*i + (shape[i],) + (1,)*(dim-i-1))
# Fill cepstrum field (all non-zero modes)
no_zero_modes = (slice(1, None),) * dim
......@@ -111,10 +107,12 @@ class AmplitudeModel(Operator):
dof_space = qht.domain[0]
sym = SymmetrizingOperator(logk_space)
phi_mean = np.array([sm, im])
phi_mean = np.array([sm, im+sm*logk_space.t_0[0]])
phi_sig = np.array([sv, iv])
self._slope = SlopeOperator(logk_space, phi_sig)
self._slope = SlopeOperator(logk_space)
self._slope = self._slope(makeOp(Field.from_global_data(
self._slope.domain,phi_sig)))
self._norm_phi_mean = Field.from_global_data(self._slope.domain,
phi_mean/phi_sig)
......
......@@ -25,6 +25,7 @@ from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import FieldAdapter
from ..operators.scaling_operator import ScalingOperator
def CorrelatedField(s_space, amplitude_model):
......@@ -42,7 +43,10 @@ def CorrelatedField(s_space, amplitude_model):
p_space = amplitude_model.target[0]
power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_model)
return ht(A*FieldAdapter(MultiDomain.make({"xi": h_space}), "xi"))
vol = h_space.scalar_dvol
#vol = 1.
vol = ScalingOperator(vol ** (-0.5),h_space)
return ht(vol(A*FieldAdapter(MultiDomain.make({"xi": h_space}), "xi")))
def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
......
......@@ -46,27 +46,29 @@ class SlopeOperator(LinearOperator):
sigmas : np.array, shape=(2,)
The slope variance and the y-intercept variance.
"""
def __init__(self, target, sigmas):
def __init__(self, target):
if not isinstance(target, LogRGSpace):
raise TypeError
self._domain = DomainTuple.make(UnstructuredDomain((2,)))
self._target = DomainTuple.make(target)
self._capability = self.TIMES | self.ADJOINT_TIMES
self._sigmas = sigmas
self.ndim = len(self.target[0].shape)
self.pos = np.zeros((self.ndim,) + self.target[0].shape)
if self.ndim == 1:
self.pos[0] = self.target[0].get_k_length_array().to_global_data()
else:
shape = self.target[0].shape
for i in range(self.ndim):
rng = np.arange(target.shape[i])
tmp = np.minimum(
rng, target.shape[i]+1-rng) * target.bindistances[i]
self.pos[i] += tmp.reshape(
(1,)*i + (shape[i],) + (1,)*(self.ndim-i-1))
if self.ndim != 1:
raise ValueError("Slope Operator only works for ndim == 1")
# Prepare pos
self.pos = np.zeros((self.ndim,) + self.target[0].shape)
shape = self.target.shape
dist = self.target[0].bindistances
for i in range(self.ndim):
ks = np.zeros(shape[i])
ks[1:] = np.minimum(shape[i] - 1 - np.arange(shape[i]-1), np.arange(shape[i]-1)) * dist[i]
self.pos[i] += ks.reshape((1,)*i + (shape[i],) + (1,)*(self.ndim-i-1))
def apply(self, x, mode):
self._check_input(x, mode)
......@@ -74,15 +76,16 @@ class SlopeOperator(LinearOperator):
# Times
if mode == self.TIMES:
inp = x.to_global_data()
res = self._sigmas[-1] * inp[-1]
res = inp[-1]
for i in range(self.ndim):
res += self._sigmas[i] * inp[i] * self.pos[i]
res = res + inp[i] * self.pos[i]
res[0] = 0.
return Field.from_global_data(self.target, res)
# Adjoint times
res = np.zeros(self.domain[0].shape, dtype=x.dtype)
xglob = x.to_global_data()
res[-1] = np.sum(xglob) * self._sigmas[-1]
res[-1] = np.sum(xglob[1:])
for i in range(self.ndim):
res[i] = np.sum(self.pos[i] * xglob) * self._sigmas[i]
res[i] = np.sum((self.pos[i] * xglob)[1:])
return Field.from_global_data(self.domain, res)
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