Commit 429b90ad authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_5' into cg_reset

parents ea3f3a14 e52b09a7
......@@ -30,7 +30,7 @@ from ..sugar import makeOp
class InverseGammaModel(Operator):
def __init__(self, domain, alpha, q):
def __init__(self, domain, alpha, q, delta=0.001):
"""Model which transforms a Gaussian into an inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows:
......@@ -42,6 +42,9 @@ class InverseGammaModel(Operator):
The mean of the pdf is at q / (alpha - 1) if alpha > 1.
The mode is q / (alpha + 1).
This transformation is implemented as a linear interpolation which
maps a Gaussian onto a inverse gamma distribution.
Parameters
----------
domain : Domain, tuple of Domain or DomainTuple
......@@ -51,30 +54,38 @@ class InverseGammaModel(Operator):
The alpha-parameter of the inverse-gamma distribution.
q : float
The q-parameter of the inverse-gamma distribution.
delta : float
distance between sampling points for linear interpolation.
"""
self._domain = self._target = DomainTuple.make(domain)
self._alpha = alpha
self._q = q
self._alpha, self._q, self._delta = float(alpha), float(q), float(delta)
self._xmin, self._xmax = -8.2, 8.2
# Precompute
xs = np.arange(self._xmin, self._xmax+2*delta, delta)
self._table = np.log(invgamma.ppf(norm.cdf(xs), self._alpha,
scale=self._q))
self._deriv = (self._table[1:]-self._table[:-1]) / delta
def apply(self, x):
self._check_input(x)
lin = isinstance(x, Linearization)
val = x.val.local_data if lin else x.local_data
# MR FIXME?!
points = np.clip(val, None, 8.2)
points = invgamma.ppf(norm.cdf(points), self._alpha, scale=self._q)
points = Field.from_local_data(self._domain, points)
val = (np.clip(val, self._xmin, self._xmax) - self._xmin) / self._delta
# Operator
fi = np.floor(val).astype(int)
w = val - fi
res = np.exp((1 - w)*self._table[fi] + w*self._table[fi + 1])
points = Field.from_local_data(self._domain, res)
if not lin:
return points
inner = norm.pdf(val)
outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(val),
self._alpha,
scale=self._q),
self._alpha, scale=self._q)
# FIXME
outer_inv = np.clip(outer_inv, 1e-20, None)
outer = 1/outer_inv
jac = makeOp(Field.from_local_data(self._domain, inner*outer))
# Derivative of linear interpolation
der = self._deriv[fi]*res
jac = makeOp(Field.from_local_data(self._domain, der))
jac = jac(x.jac)
return x.new(points, jac)
......
......@@ -75,14 +75,13 @@ class FFTOperator(LinearOperator):
target.check_codomain(adom)
def apply(self, x, mode):
from pyfftw.interfaces.numpy_fft import fftn, ifftn
self._check_input(x, mode)
ncells = x.domain[self._space].size
if x.domain[self._space].harmonic: # harmonic -> position
func = fftn
func = fft.fftn
fct = 1.
else:
func = ifftn
func = fft.ifftn
fct = ncells
axes = x.domain.axes[self._space]
tdom = self._tgt(mode)
......
......@@ -219,7 +219,7 @@ class _OpSum(Operator):
lin1 = self._op1(Linearization.make_var(v1, wm))
lin2 = self._op2(Linearization.make_var(v2, wm))
op = lin1._jac._myadd(lin2._jac, False)
res = lin1.new(lin1._val+lin2._val, op(x.jac))
res = lin1.new(lin1._val.unite(lin2._val), op(x.jac))
if lin1._metric is not None and lin2._metric is not None:
res = res.add_metric(lin1._metric + lin2._metric)
return 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