Commit 0e63c553 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_5' into clipped_exp

parents eb80cb32 e52b09a7
......@@ -24,6 +24,7 @@ from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import DOFDistributor, PowerDistributor
from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
from .operators.contraction_operator import ContractionOperator
from .operators.linear_interpolation import LinearInterpolator
from .operators.endomorphic_operator import EndomorphicOperator
from .operators.exp_transform import ExpTransform
from .operators.harmonic_operators import (
......
......@@ -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)
......
from __future__ import absolute_import, division, print_function
import numpy as np
......@@ -101,6 +102,12 @@ class Linearization(object):
def __rtruediv__(self, other):
return self.inverse().__mul__(other)
def __pow__(self, power):
if not np.isscalar(power):
return NotImplemented
return self.new(self._val**power,
makeOp(self._val**(power-1)).scale(power)(self._jac))
def inverse(self):
return self.new(1./self._val, makeOp(-1./(self._val**2))(self._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)
......
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from __future__ import absolute_import, division, print_function
from ..compat import *
from .. import Field, UnstructuredDomain
from ..sugar import makeDomain
from .linear_operator import LinearOperator
from numpy import (array, prod, mgrid, int64, arange, ravel_multi_index, zeros,
abs, ravel)
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import aslinearoperator
class LinearInterpolator(LinearOperator):
def __init__(self, domain, positions):
"""
:param domain:
RGSpace
:param target:
UnstructuredDomain, shape (ndata,)
:param positions:
positions at which to interpolate
Field with UnstructuredDomain, shape (dim, ndata)
"""
self._domain = makeDomain(domain)
N_points = positions.shape[1]
self._target = makeDomain(UnstructuredDomain(N_points))
self._capability = self.TIMES | self.ADJOINT_TIMES
self._build_mat(positions, N_points)
def _build_mat(self, positions, N_points):
ndim = positions.shape[0]
mg = mgrid[(slice(0, 2),)*ndim]
mg = array(list(map(ravel, mg)))
dist = array(self.domain[0].distances).reshape((-1, 1))
pos = positions/dist
excess = pos-pos.astype(int64)
pos = pos.astype(int64)
data = zeros((len(mg[0]), N_points))
ii = zeros((len(mg[0]), N_points), dtype=int64)
jj = zeros((len(mg[0]), N_points), dtype=int64)
for i in range(len(mg[0])):
factor = prod(abs(1-mg[:, i].reshape((-1, 1))-excess), axis=0)
data[i, :] = factor
fromi = pos+mg[:, i].reshape((-1, 1))
ii[i, :] = arange(N_points)
jj[i, :] = ravel_multi_index(fromi, self.domain.shape)
self._mat = coo_matrix((data.reshape(-1),
(ii.reshape(-1), jj.reshape(-1))),
(N_points, prod(self.domain.shape)))
self._mat = aslinearoperator(self._mat)
def apply(self, x, mode):
self._check_input(x, mode)
x_val = x.to_global_data()
if mode == self.TIMES:
res = self._mat.matvec(x_val.reshape((-1,)))
return Field.from_global_data(self.target, res)
res = self._mat.rmatvec(x_val).reshape(self.domain.shape)
return Field.from_global_data(self.domain, res)
# import numpy as np
# from ..domains.rg_space import RGSpace
# import itertools
#
#
# class LinearInterpolator(LinearOperator):
# def __init__(self, domain, positions):
# """
# :param domain:
# RGSpace
# :param target:
# UnstructuredDomain, shape (ndata,)
# :param positions:
# positions at which to interpolate
# Field with UnstructuredDomain, shape (dim, ndata)
# """
# if not isinstance(domain, RGSpace):
# raise TypeError("RGSpace needed")
# if np.any(domain.shape < 2):
# raise ValueError("RGSpace shape too small")
# if positions.ndim != 2:
# raise ValueError("positions must be a 2D array")
# ndim = len(domain.shape)
# if positions.shape[0] != ndim:
# raise ValueError("shape mismatch")
# self._domain = makeDomain(domain)
# N_points = positions.shape[1]
# dist = np.array(domain.distances).reshape((ndim, -1))
# self._pos = positions/dist
# shp = np.array(domain.shape, dtype=int).reshape((ndim, -1))
# self._idx = np.maximum(0, np.minimum(shp-2, self._pos.astype(int)))
# self._pos -= self._idx
# tmp = tuple([0, 1] for i in range(ndim))
# self._corners = np.array(list(itertools.product(*tmp)))
# self._target = makeDomain(UnstructuredDomain(N_points))
# self._capability = self.TIMES | self.ADJOINT_TIMES
#
# def apply(self, x, mode):
# self._check_input(x, mode)
# x = x.to_global_data()
# ndim = len(self._domain.shape)
#
# res = np.zeros(self._tgt(mode).shape, dtype=x.dtype)
# for corner in self._corners:
# corner = corner.reshape(ndim, -1)
# idx = self._idx+corner
# idx2 = tuple(idx[t, :] for t in range(idx.shape[0]))
# wgt = np.prod(self._pos*corner+(1-self._pos)*(1-corner), axis=0)
# if mode == self.TIMES:
# res += wgt*x[idx2]
# else:
# np.add.at(res, idx2, wgt*x)
# return Field.from_global_data(self._tgt(mode), res)
from __future__ import absolute_import, division, print_function
import numpy as np
from ..compat import *
from ..utilities import NiftyMetaBase
......@@ -29,7 +30,7 @@ class Operator(NiftyMetaBase()):
s = "The operator's and field's domains don't match."
from ..domain_tuple import DomainTuple
from ..multi_domain import MultiDomain
if not isinstance(dom_op, [DomainTuple, MultiDomain]):
if not isinstance(dom_op, (DomainTuple, MultiDomain,)):
s += " Your operator's domain is neither a `DomainTuple`" \
" nor a `MultiDomain`."
raise ValueError(s)
......@@ -67,6 +68,16 @@ class Operator(NiftyMetaBase()):
return NotImplemented
return _OpSum(self, x)
def __sub__(self, x):
if not isinstance(x, Operator):
return NotImplemented
return _OpSum(self, -x)
def __pow__(self, power):
if not np.isscalar(power):
return NotImplemented
return _OpChain.make((_PowerOp(self.target, power), self))
def apply(self, x):
raise NotImplementedError
......@@ -108,7 +119,17 @@ class _FunctionApplier(Operator):
return getattr(x, self._funcname)()
# FIXME Is this class used except in _OpChain? Can it be merged?
class _PowerOp(Operator):
def __init__(self, domain, power):
from ..sugar import makeDomain
self._domain = self._target = makeDomain(domain)
self._power = power
def apply(self, x):
self._check_input(x)
return x**self._power
class _CombinedOperator(Operator):
def __init__(self, ops, _callingfrommake=False):
if not _callingfrommake:
......@@ -198,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
......@@ -58,12 +58,13 @@ class Consistency_Tests(unittest.TestCase):
op = a+b
ift.extra.consistency_check(op, dtype, dtype)
@expand(product(_h_spaces + _p_spaces + _pow_spaces,
[np.float64, np.complex128]))
def testVdotOperator(self, sp, dtype):
op = ift.VdotOperator(ift.Field.from_random("normal", sp,
dtype=dtype))
ift.extra.consistency_check(op, dtype, dtype)
def testLinearInterpolator(self):
sp = ift.RGSpace((10,8), distances=(0.1, 3.5))
pos = np.random.rand(2, 23)
pos[0,:] *= 0.9
pos[1,:] *= 7*3.5
op = ift.LinearInterpolator(sp, pos)
ift.extra.consistency_check(op)
@expand(product([(ift.RGSpace(10, harmonic=True), 4, 0),
(ift.RGSpace((24, 31), distances=(0.4, 2.34),
......
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