Commit 50168dda authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'UniformOperator' into 'NIFTy_6'

Uniform operator

See merge request !419
parents 3541b181 c7b89f8a
Pipeline #70801 passed with stages
in 29 minutes and 25 seconds
......@@ -72,7 +72,7 @@ from .sugar import *
from .plot import Plot
from .library.special_distributions import InverseGammaOperator
from .library.special_distributions import InverseGammaOperator, UniformOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator,
dynamic_lightcone_operator)
......
......@@ -17,44 +17,82 @@
import numpy as np
from scipy.stats import invgamma, norm
from scipy.interpolate import CubicSpline
from .. import Adder
from ..domain_tuple import DomainTuple
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..linearization import Linearization
from ..operators.operator import Operator
from ..sugar import makeOp
def _f_on_np(f, arr):
fld = Field.from_raw(UnstructuredDomain(arr.shape), arr)
return f(fld).val
class _InterpolationOperator(Operator):
def __init__(self, domain, func, xmin, xmax, delta, table_func=None, inverse_table_func=None):
"""
Calculates a function pointwise on a field by interpolation.
Can be supplied with a table_func to increase the interpolation accuracy,
Best results are achieved when table_func(func) is roughly linear.
Parameters
----------
domain : Domain, tuple of Domain or DomainTuple
The domain on which the field shall be defined. This is at the same
time the domain and the target of the operator.
func : function
The function which is applied on the field. Assumed to act on numpy
arrays.
xmin : float
The smallest value for which func will be evaluated.
xmax : float
The largest value for which func will be evaluated.
delta : float
Distance between sampling points for linear interpolation.
table_func : function
Non-linear function applied to table in order to transform the table
to a more linear space. Assumed to act on `Linearization`s, optional.
inv_table_func : function
Inverse of table_func, optional.
"""
def __init__(self, domain, func, xmin, xmax, delta, table_func=None, inv_table_func=None):
self._domain = self._target = DomainTuple.make(domain)
self._xmin, self._xmax = float(xmin), float(xmax)
self._d = float(delta)
self._xs = np.arange(xmin, xmax+2*self._d, self._d)
self._xs = np.arange(xmin, xmax, self._d)
self._table = func(self._xs)
if table_func is not None:
foo = func(np.random.randn(10))
np.testing.assert_allclose(foo, inverse_table_func(table_func(foo)))
self._table = table_func(func(self._xs))
self._deriv = (self._table[1:]-self._table[:-1]) / self._d
self._inv_table_func = inverse_table_func
if inv_table_func is None:
raise ValueError
a = func(np.random.randn(10))
a1 = _f_on_np(lambda x: inv_table_func(table_func(x)), a)
np.testing.assert_allclose(a, a1)
self._table = _f_on_np(table_func, self._table)
self._interpolator = CubicSpline(self._xs, self._table)
self._deriv = self._interpolator.derivative()
self._inv_table_func = inv_table_func
def apply(self, x):
self._check_input(x)
lin = isinstance(x, Linearization)
xval = x.val.val if lin else x.val
val = (np.clip(xval, self._xmin, self._xmax) - self._xmin) / self._d
fi = np.floor(val).astype(int)
w = val - fi
res = self._inv_table_func((1-w)*self._table[fi] + w*self._table[fi+1])
resfld = Field(self._domain, res)
if not lin:
return resfld
jac = makeOp(Field(self._domain, self._deriv[fi]*res))
return x.new(resfld, jac)
res = self._interpolator(xval)
res = Field(self._domain, res)
if lin:
res = x.new(res, makeOp(Field(self._domain, self._deriv(xval))))
if self._inv_table_func is not None:
res = self._inv_table_func(res)
return res
def InverseGammaOperator(domain, alpha, q, delta=0.001):
"""Transforms a Gaussian into an inverse gamma distribution.
def InverseGammaOperator(domain, alpha, q, delta=1e-2):
"""Transforms a Gaussian with unit covariance and zero mean into an
inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows:
......@@ -67,7 +105,7 @@ def InverseGammaOperator(domain, alpha, q, delta=0.001):
The mode is :math:`q / (\\alpha + 1)`.
This transformation is implemented as a linear interpolation which maps a
Gaussian onto a inverse gamma distribution.
Gaussian onto an inverse gamma distribution.
Parameters
----------
......@@ -76,14 +114,57 @@ def InverseGammaOperator(domain, alpha, q, delta=0.001):
time the domain and the target of the operator.
alpha : float
The alpha-parameter of the inverse-gamma distribution.
q : float
q : float or Field
The q-parameter of the inverse-gamma distribution.
delta : float
distance between sampling points for linear interpolation.
Distance between sampling points for linear interpolation.
"""
op = _InterpolationOperator(domain,
lambda x: invgamma.ppf(norm.cdf(x), float(alpha)),
-8.2, 8.2, delta, np.log, np.exp)
op = _InterpolationOperator(domain, lambda x: invgamma.ppf(norm._cdf(x), float(alpha)),
-8.2, 8.2, delta, lambda x: x.log(), lambda x: x.exp())
if np.isscalar(q):
return op.scale(q)
return makeOp(q) @ op
class UniformOperator(Operator):
"""
Transforms a Gaussian with unit covariance and zero mean into a uniform
distribution. The uniform distribution's support is ``[loc, loc + scale]``.
Parameters
----------
domain : Domain, tuple of Domain or DomainTuple
The domain on which the field shall be defined. This is at the same
time the domain and the target of the operator.
loc: float
scale: float
"""
def __init__(self, domain, loc=0, scale=1):
self._target = self._domain = DomainTuple.make(domain)
self.loc = loc
self.scale = scale
def apply(self, x):
self._check_input(x)
lin = isinstance(x, Linearization)
xval = x.val.val if lin else x.val
res = Field(self._target, norm._cdf(xval) * self.scale + self.loc)
if not lin:
return res
jac = makeOp(Field(self._domain, norm._pdf(xval)*self.scale))
return x.new(res, jac)
@staticmethod
def uni(field, loc=0, high=1):
foo = norm._cdf(field.val) * scale + loc
return Field(field.domain, foo)
@staticmethod
def inverse_uni(field, low=0, scale=1):
res = norm._ppf(field.val/scale - loc)
return Field(field.domain, res)
......@@ -32,7 +32,7 @@ space = list2fixture([ift.GLSpace(15),
seed = list2fixture([4, 78, 23])
def testInverseGammaAccuracy(space, seed):
def testInterpolationAccuracy(space, seed):
S = ift.ScalingOperator(space, 1.)
pos = S.draw_sample()
alpha = 1.5
......@@ -45,3 +45,4 @@ def testInverseGammaAccuracy(space, seed):
arr1 = op(pos).val
arr0 = invgamma.ppf(norm.cdf(pos.val), alpha, scale=q)
assert_allclose(arr0, arr1)
......@@ -89,15 +89,16 @@ def testBinary(type1, type2, space, seed):
ift.extra.check_jacobian_consistency(model, pos, ntries=20)
def testPointModel(space, seed):
def testSpecialDistributionOps(space, seed):
np.random.seed(seed)
S = ift.ScalingOperator(space, 1.)
pos = S.draw_sample()
alpha = 1.5
q = 0.73
model = ift.InverseGammaOperator(space, alpha, q)
# FIXME All those cdfs and ppfs are not very accurate
ift.extra.check_jacobian_consistency(model, pos, tol=1e-2, ntries=20)
ift.extra.check_jacobian_consistency(model, pos, ntries=20)
model = ift.UniformOperator(space, alpha, q)
ift.extra.check_jacobian_consistency(model, pos, ntries=20)
@pmp('neg', [True, False])
......
Supports Markdown
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