Commit 97829024 authored by Martin Reinecke's avatar Martin Reinecke

tweaks; added some FIXMEs

parent 966f258b
......@@ -104,9 +104,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
for i in range(dim):
ks = np.minimum(shape[i] - np.arange(shape[i]) +
1, np.arange(shape[i])) * dist[i]
fst_dims = (1,) * i
lst_dims = (1,) * (dim - i - 1)
q_array[i] += ks.reshape(fst_dims + (shape[i],) + lst_dims)
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
......@@ -117,10 +115,9 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
# Fill cepstrum field (zero-mode subspaces)
for i in range(dim):
# Prepare indices
fst_dims = (slice(None),) * i
lst_dims = (slice(None),) * (dim - i - 1)
sl = fst_dims + (slice(1, None),) + lst_dims
sl2 = fst_dims + (0,) + lst_dims
fst_dims = (slice(None),)*i
sl = fst_dims + (slice(1, None),)
sl2 = fst_dims + (0,)
# Do summation
cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
......
......@@ -35,10 +35,10 @@ class BernoulliEnergy(Energy):
self._d = d
p_val = self._p.value
self._value = -self._d.vdot(log(p_val)) - (1. - d).vdot(log(1.-p_val))
self._value = -self._d.vdot(log(p_val)) - (1.-d).vdot(log(1.-p_val))
if isnan(self._value):
self._value = inf
metric = makeOp(1./((p_val) * (1.-p_val)))
metric = makeOp(1. / (p_val * (1.-p_val)))
self._gradient = self._p.jacobian.adjoint_times(metric(p_val-d))
self._metric = SandwichOperator.make(self._p.jacobian, metric)
......
from ..operators.fft_operator import FFTOperator
from ..field import Field
from ..multi.multi_field import MultiField
from ..models.local_nonlinearity import PointwiseExponential
from ..operators.power_distributor import PowerDistributor
from ..models.variable import Variable
from ..domain_tuple import DomainTuple
from ..operators.domain_distributor import DomainDistributor
from ..operators.harmonic_transform_operator \
import HarmonicTransformOperator
def make_correlated_field(s_space, amplitude_model):
'''
Method for construction of correlated fields
......@@ -8,21 +19,14 @@ def make_correlated_field(s_space, amplitude_model):
amplitude_model : model for correlation structure
'''
from ..operators.fft_operator import FFTOperator
from ..field import Field
from ..multi.multi_field import MultiField
from ..models.local_nonlinearity import PointwiseExponential
from ..operators.power_distributor import PowerDistributor
from ..models.variable import Variable
h_space = s_space.get_default_codomain()
ht = FFTOperator(h_space, s_space)
p_space = amplitude_model.value.domain[0]
power_distributor = PowerDistributor(h_space, p_space)
position = {}
position['xi'] = Field.from_random('normal', h_space)
position['tau'] = amplitude_model.position['tau']
position['phi'] = amplitude_model.position['phi']
position = MultiField.from_dict(position)
position = MultiField.from_dict({
'xi': Field.from_random('normal', h_space),
'tau': amplitude_model.position['tau'],
'phi': amplitude_model.position['phi']})
xi = Variable(position)['xi']
A = power_distributor(amplitude_model)
......@@ -39,16 +43,6 @@ def make_mf_correlated_field(s_space_spatial, s_space_energy,
'''
Method for construction of correlated multi-frequency fields
'''
from ..operators.fft_operator import FFTOperator
from ..field import Field
from ..multi.multi_field import MultiField
from ..models.local_nonlinearity import PointwiseExponential
from ..operators.power_distributor import PowerDistributor
from ..models.variable import Variable
from ..domain_tuple import DomainTuple
from ..operators.domain_distributor import DomainDistributor
from ..operators.harmonic_transform_operator \
import HarmonicTransformOperator
h_space_spatial = s_space_spatial.get_default_codomain()
h_space_energy = s_space_energy.get_default_codomain()
h_space = DomainTuple.make((h_space_spatial, h_space_energy))
......
......@@ -21,6 +21,7 @@ from ..operators.sandwich_operator import SandwichOperator
from ..utilities import memo
# MR FIXME documentation incomplete
class GaussianEnergy(Energy):
def __init__(self, inp, mean=None, covariance=None):
"""
......@@ -39,30 +40,29 @@ class GaussianEnergy(Energy):
@property
@memo
def residual(self):
if self._mean is not None:
return self._inp.value - self._mean
return self._inp.value
def _residual(self):
return self._inp.value if self._mean is None else \
self._inp.value - self._mean
@property
@memo
def _icovres(self):
return self._residual if self._cov is None else \
self._cov.inverse_times(self._residual)
@property
@memo
def value(self):
if self._cov is None:
return .5 * self.residual.vdot(self.residual).real
return .5 * self.residual.vdot(
self._cov.inverse_times(self.residual)).real
return .5 * self._residual.vdot(self._icovres).real
@property
@memo
def gradient(self):
if self._cov is None:
return self._inp.jacobian.adjoint_times(self.residual)
return self._inp.jacobian.adjoint_times(
self._cov.inverse_times(self.residual))
return self._inp.jacobian.adjoint_times(self._icovres)
@property
@memo
def metric(self):
if self._cov is None:
return SandwichOperator.make(self._inp.jacobian, None)
return SandwichOperator.make(self._inp.jacobian, self._cov.inverse)
return SandwichOperator.make(
self._inp.jacobian,
None if self._cov is None else self._cov.inverse)
......@@ -50,6 +50,7 @@ class PointSources(Model):
foo = invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q)
return Field.from_local_data(field.domain, foo)
# MR FIXME: is this function needed?
@staticmethod
def IG_prime(field, alpha, q):
inner = norm.pdf(field.local_data)
......@@ -57,9 +58,9 @@ class PointSources(Model):
scale=q), alpha, scale=q)
# # FIXME
# outer = np.clip(outer, 1e-20, None)
outer = 1/outer
return Field.from_local_data(field.domain, inner*outer)
return Field.from_local_data(field.domain, inner/outer)
# MR FIXME: why does this take an np.ndarray instead of a Field?
@staticmethod
def inverseIG(u, alpha, q):
res = norm.ppf(invgamma.cdf(u, alpha, scale=q))
......
......@@ -42,7 +42,6 @@ class PoissonianEnergy(Energy):
self._value = inf
self._gradient = self._lamb.jacobian.adjoint_times(1 - d/lamb_val)
# metric = makeOp(d/lamb_val/lamb_val)
metric = makeOp(1./lamb_val)
self._metric = SandwichOperator.make(self._lamb.jacobian, metric)
......
......@@ -44,7 +44,8 @@ class BlockDiagonalOperator(EndomorphicOperator):
def _combine_chain(self, op):
if self._domain is not op._domain:
raise ValueError("domain mismatch")
res = {key : v1*v2 for key, v1, v2 in zip(self._domain.keys(), self._ops, op._ops)}
res = {key: v1*v2
for key, v1, v2 in zip(self._domain.keys(), self._ops, op._ops)}
return BlockDiagonalOperator(self._domain, res)
def _combine_sum(self, op, selfneg, opneg):
......
......@@ -66,6 +66,6 @@ class MultiDomain(object):
def __str__(self):
res = "MultiDomain:\n"
for key, dom in zip(self._keys, self._domains):
for key, dom in zip(self._keys, self._domains):
res += key+": "+str(dom)+"\n"
return res
......@@ -30,4 +30,4 @@ def MultiAdaptor(target):
"""
if not isinstance(target, MultiDomain) or len(target) > 1:
raise TypeError
return SelectionOperator(target,target.keys()[0]).adjoint
return SelectionOperator(target, target.keys()[0]).adjoint
......@@ -233,8 +233,8 @@ def makeOp(input):
if isinstance(input, Field):
return DiagonalOperator(input)
if isinstance(input, MultiField):
return BlockDiagonalOperator(input.domain, {key: makeOp(val)
for key, val in input.items()})
return BlockDiagonalOperator(
input.domain, {key: makeOp(val) for key, val in input.items()})
raise NotImplementedError
# Arithmetic functions working on Fields
......
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