Commit 35cf799b authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'new_apply' into 'NIFTy_6'

New apply

See merge request !418
parents a50ec021 646ccaa7
Pipeline #70644 passed with stages
in 17 minutes and 51 seconds
Changes since NIFTy 5:
Interface Change for non-linear Operators
=========================================
The method `Operator.apply()` takes a `Linearization` or a `Field` or a
`MultiField` as input. This has not changed. However, now each non-linear
operator assumes that the input `Linearization` comes with an identity operator
as jacobian. Also it is assumed that the `apply()` method returns a
`Linearization` with the jacobian of the operator itself. The user is not in
charge anymore of stacking together the jacobians of operator chains. Something
like `x.jac` should not appear in any self-written `apply()` methods. The method
`Operator._check_input` tests if this condition is met. The same goes for the
metric. There is no need anymore to call `SandwichOperator` in an `apply()`
method when implementing new energies. This change should not lead to unexpected
behaviour since both `Operator._check_input()` and
`extra.check_jacobian_consistency()` tests for the new conditions to be
fulfilled.
Special functions for complete Field reduction operations
=========================================================
So far, reduction operations called on Fields (like `vdot`, `sum`, `integrate`,
`mean`, `var`, `std`, `prod` etc.) returned a scalar when the reduction was
carried out over all domains, and otherwise a `Field`.
Having the data type of the returned value depend on input parameters is
extremely confusing, so all of these reduction operations now always return a
Field. We also introduced another set of reduction operations which always
operate over all subdomains and therefore don't take a `spaces` argument; they
are named `s_vdot`, `s_sum` etc. and always return a scalar.
Updates regarding correlated fields
===================================
......
......@@ -67,7 +67,7 @@ if __name__ == '__main__':
# Compute likelihood and Hamiltonian
position = ift.from_random('normal', harmonic_space)
likelihood = ift.BernoulliEnergy(data)(p)
likelihood = ift.BernoulliEnergy(data) @ p
ic_newton = ift.DeltaEnergyController(
name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
minimizer = ift.NewtonCG(ic_newton)
......
......@@ -96,7 +96,7 @@ if __name__ == '__main__':
data = lamb(mock_position)
data = np.random.poisson(data.val.astype(np.float64))
data = ift.Field.from_raw(d_space, data)
likelihood = ift.PoissonianEnergy(data)(lamb)
likelihood = ift.PoissonianEnergy(data) @ lamb
# Settings for minimization
ic_newton = ift.DeltaEnergyController(
......
......@@ -122,8 +122,7 @@ if __name__ == '__main__':
N_samples = 20
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(
mean=data, inverse_covariance=N.inverse)(signal_response)
likelihood = ift.GaussianEnergy(mean=data, inverse_covariance=N.inverse) @ signal_response
H = ift.StandardHamiltonian(likelihood, ic_sampling)
# Begin minimization
......
......@@ -81,63 +81,63 @@ class PolynomialResponse(ift.LinearOperator):
return ift.makeField(self._tgt(mode), out)
# Generate some mock data
N_params = 10
N_samples = 100
size = (12,)
x = np.random.random(size) * 10
y = np.sin(x**2) * x**3
var = np.full_like(y, y.var() / 10)
var[-2] *= 4
var[5] /= 2
y[5] -= 0
# Set up minimization problem
p_space = ift.UnstructuredDomain(N_params)
params = ift.full(p_space, 0.)
R = PolynomialResponse(p_space, x)
ift.extra.consistency_check(R)
d_space = R.target
d = ift.makeField(d_space, y)
N = ift.DiagonalOperator(ift.makeField(d_space, var))
IC = ift.DeltaEnergyController(tol_rel_deltaE=1e-12, iteration_limit=200)
likelihood = ift.GaussianEnergy(d, N)(R)
Ham = ift.StandardHamiltonian(likelihood, IC)
H = ift.EnergyAdapter(params, Ham, want_metric=True)
# Minimize
minimizer = ift.NewtonCG(IC)
H, _ = minimizer(H)
# Draw posterior samples
metric = Ham(ift.Linearization.make_var(H.position, want_metric=True)).metric
samples = [metric.draw_sample(from_inverse=True) + H.position
for _ in range(N_samples)]
# Plotting
plt.errorbar(x, y, np.sqrt(var), fmt='ko', label='Data with error bars')
xmin, xmax = x.min(), x.max()
xs = np.linspace(xmin, xmax, 100)
sc = ift.StatCalculator()
for ii in range(len(samples)):
sc.add(samples[ii])
ys = polynomial(samples[ii], xs)
if ii == 0:
plt.plot(xs, ys, 'k', alpha=.05, label='Posterior samples')
continue
plt.plot(xs, ys, 'k', alpha=.05)
ys = polynomial(H.position, xs)
plt.plot(xs, ys, 'r', linewidth=2., label='Interpolation')
plt.legend()
plt.savefig('fit.png')
plt.close()
# Print parameters
mean = sc.mean.val
sigma = np.sqrt(sc.var.val)
for ii in range(len(mean)):
print('Coefficient x**{}: {:.2E} +/- {:.2E}'.format(ii, mean[ii],
sigma[ii]))
if __name__ == '__main__':
# Generate some mock data
N_params = 10
N_samples = 100
size = (12,)
x = np.random.random(size) * 10
y = np.sin(x**2) * x**3
var = np.full_like(y, y.var() / 10)
var[-2] *= 4
var[5] /= 2
y[5] -= 0
# Set up minimization problem
p_space = ift.UnstructuredDomain(N_params)
params = ift.full(p_space, 0.)
R = PolynomialResponse(p_space, x)
ift.extra.consistency_check(R)
d_space = R.target
d = ift.makeField(d_space, y)
N = ift.DiagonalOperator(ift.makeField(d_space, var))
IC = ift.DeltaEnergyController(tol_rel_deltaE=1e-12, iteration_limit=200)
likelihood = ift.GaussianEnergy(d, N) @ R
Ham = ift.StandardHamiltonian(likelihood, IC)
H = ift.EnergyAdapter(params, Ham, want_metric=True)
# Minimize
minimizer = ift.NewtonCG(IC)
H, _ = minimizer(H)
# Draw posterior samples
metric = Ham(ift.Linearization.make_var(H.position, want_metric=True)).metric
samples = [metric.draw_sample(from_inverse=True) + H.position
for _ in range(N_samples)]
# Plotting
plt.errorbar(x, y, np.sqrt(var), fmt='ko', label='Data with error bars')
xmin, xmax = x.min(), x.max()
xs = np.linspace(xmin, xmax, 100)
sc = ift.StatCalculator()
for ii in range(len(samples)):
sc.add(samples[ii])
ys = polynomial(samples[ii], xs)
if ii == 0:
plt.plot(xs, ys, 'k', alpha=.05, label='Posterior samples')
continue
plt.plot(xs, ys, 'k', alpha=.05)
ys = polynomial(H.position, xs)
plt.plot(xs, ys, 'r', linewidth=2., label='Interpolation')
plt.legend()
plt.savefig('fit.png')
plt.close()
# Print parameters
mean = sc.mean.val
sigma = np.sqrt(sc.var.val)
for ii in range(len(mean)):
print('Coefficient x**{}: {:.2E} +/- {:.2E}'.format(ii, mean[ii], sigma[ii]))
......@@ -15,7 +15,10 @@ from .multi_domain import MultiDomain
from .field import Field
from .multi_field import MultiField
from .operators.operator import Operator
from .operators.linear_operator import LinearOperator
from .operators.adder import Adder
from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import DOFDistributor, PowerDistributor
......@@ -28,7 +31,6 @@ from .operators.harmonic_operators import (
HarmonicSmoothingOperator)
from .operators.field_zero_padder import FieldZeroPadder
from .operators.inversion_enabler import InversionEnabler
from .operators.linear_operator import LinearOperator
from .operators.mask_operator import MaskOperator
from .operators.regridding_operator import RegriddingOperator
from .operators.sampling_enabler import SamplingEnabler
......@@ -67,6 +69,7 @@ from .minimization.energy_adapter import EnergyAdapter
from .minimization.metric_gaussian_kl import MetricGaussianKL
from .sugar import *
from .plot import Plot
from .library.special_distributions import InverseGammaOperator
......
......@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -125,7 +125,7 @@ class LMSpace(StructuredDomain):
# evaluate the kernel function at the required thetas
kernel_sphere = Field.from_raw(gl, func(theta))
# normalize the kernel such that the integral over the sphere is 4pi
kernel_sphere = kernel_sphere * (4 * np.pi / kernel_sphere.integrate())
kernel_sphere = kernel_sphere * (4 * np.pi / kernel_sphere.s_integrate())
# compute the spherical harmonic coefficients of the kernel
op = HarmonicTransformOperator(lm0, gl)
kernel_lm = op.adjoint_times(kernel_sphere.weight(1)).val
......
......@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -168,7 +168,7 @@ class RGSpace(StructuredDomain):
op = HarmonicTransformOperator(self, self.get_default_codomain())
dist = op.target[0]._get_dist_array()
kernel = Field(op.target, func(dist.val))
kernel = kernel / kernel.integrate()
kernel = kernel / kernel.s_integrate()
return op.adjoint_times(kernel.weight(1))
def get_default_codomain(self):
......
......@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -44,8 +44,8 @@ def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol,
return
f1 = from_random("normal", op.domain, dtype=domain_dtype)
f2 = from_random("normal", op.target, dtype=target_dtype)
res1 = f1.vdot(op.adjoint_times(f2))
res2 = op.times(f1).vdot(f2)
res1 = f1.s_vdot(op.adjoint_times(f2))
res2 = op.times(f1).s_vdot(f2)
if only_r_linear:
res1, res2 = res1.real, res2.real
np.testing.assert_allclose(res1, res2, atol=atol, rtol=rtol)
......@@ -83,7 +83,7 @@ def _check_linearity(op, domain_dtype, atol, rtol):
assert_allclose(val1, val2, atol=atol, rtol=rtol)
def _actual_domain_check(op, domain_dtype=None, inp=None):
def _actual_domain_check_linear(op, domain_dtype=None, inp=None):
needed_cap = op.TIMES
if (op.capability & needed_cap) != needed_cap:
return
......@@ -98,21 +98,25 @@ def _actual_domain_check(op, domain_dtype=None, inp=None):
def _actual_domain_check_nonlinear(op, loc):
assert isinstance(loc, (Field, MultiField))
assert_(loc.domain is op.domain)
lin = Linearization.make_var(loc, False)
reslin = op(lin)
assert_(lin.domain is op.domain)
assert_(lin.target is op.domain)
assert_(lin.val.domain is lin.domain)
for wm in [False, True]:
lin = Linearization.make_var(loc, wm)
reslin = op(lin)
assert_(lin.domain is op.domain)
assert_(lin.target is op.domain)
assert_(lin.val.domain is lin.domain)
assert_(reslin.domain is op.domain)
assert_(reslin.target is op.target)
assert_(reslin.val.domain is reslin.target)
assert_(reslin.domain is op.domain)
assert_(reslin.target is op.target)
assert_(reslin.val.domain is reslin.target)
assert_(reslin.target is op.target)
assert_(reslin.jac.domain is reslin.domain)
assert_(reslin.jac.target is reslin.target)
_actual_domain_check(reslin.jac, inp=loc)
_actual_domain_check(reslin.jac.adjoint, inp=reslin.jac(loc))
assert_(reslin.target is op.target)
assert_(reslin.jac.domain is reslin.domain)
assert_(reslin.jac.target is reslin.target)
_actual_domain_check_linear(reslin.jac, inp=loc)
_actual_domain_check_linear(reslin.jac.adjoint, inp=reslin.jac(loc))
if reslin.metric is not None:
assert_(reslin.metric.domain is reslin.metric.target)
assert_(reslin.metric.domain is op.domain)
def _domain_check(op):
......@@ -140,16 +144,16 @@ def _performance_check(op, pos, raise_on_fail):
return self._count
for wm in [False, True]:
cop = CountingOp(op.domain)
op = op @ cop
op(pos)
myop = op @ cop
myop(pos)
cond = [cop.count != 1]
lin = op(2*Linearization.make_var(pos, wm))
lin = myop(2*Linearization.make_var(pos, wm))
cond.append(cop.count != 2)
lin.jac(pos)
cond.append(cop.count != 3)
lin.jac.adjoint(lin.val)
cond.append(cop.count != 4)
if wm and op.target is DomainTuple.scalar_domain():
if lin.metric is not None:
lin.metric(pos)
cond.append(cop.count != 6)
if any(cond):
......@@ -195,10 +199,10 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
if not isinstance(op, LinearOperator):
raise TypeError('This test tests only linear operators.')
_domain_check(op)
_actual_domain_check(op, domain_dtype)
_actual_domain_check(op.adjoint, target_dtype)
_actual_domain_check(op.inverse, target_dtype)
_actual_domain_check(op.adjoint.inverse, domain_dtype)
_actual_domain_check_linear(op, domain_dtype)
_actual_domain_check_linear(op.adjoint, target_dtype)
_actual_domain_check_linear(op.inverse, target_dtype)
_actual_domain_check_linear(op.adjoint.inverse, domain_dtype)
_check_linearity(op, domain_dtype, atol, rtol)
_check_linearity(op.adjoint, target_dtype, atol, rtol)
_check_linearity(op.inverse, target_dtype, atol, rtol)
......@@ -214,7 +218,7 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
def _get_acceptable_location(op, loc, lin):
if not np.isfinite(lin.val.sum()):
if not np.isfinite(lin.val.s_sum()):
raise ValueError('Initial value must be finite')
dir = from_random("normal", loc.domain)
dirder = lin.jac(dir)
......@@ -227,7 +231,7 @@ def _get_acceptable_location(op, loc, lin):
try:
loc2 = loc+dir
lin2 = op(Linearization.make_var(loc2, lin.want_metric))
if np.isfinite(lin2.val.sum()) and abs(lin2.val.sum()) < 1e20:
if np.isfinite(lin2.val.s_sum()) and abs(lin2.val.s_sum()) < 1e20:
break
except FloatingPointError:
pass
......@@ -281,7 +285,7 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100, perf_check=True):
dirder = linmid.jac(dir)
numgrad = (lin2.val-lin.val)
xtol = tol * dirder.norm() / np.sqrt(dirder.size)
if (abs(numgrad-dirder) <= xtol).all():
if (abs(numgrad-dirder) <= xtol).s_all():
break
dir = dir*0.5
dirnorm *= 0.5
......
......@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -284,7 +284,7 @@ class Field(object):
from .operators.outer_product_operator import OuterProduct
return OuterProduct(self, x.domain)(x)
def vdot(self, x=None, spaces=None):
def vdot(self, x, spaces=None):
"""Computes the dot product of 'self' with x.
Parameters
......@@ -312,11 +312,33 @@ class Field(object):
spaces = utilities.parse_spaces(spaces, ndom)
if len(spaces) == ndom:
return np.vdot(self._val, x._val)
return Field.scalar(np.array(np.vdot(self._val, x._val)))
# If we arrive here, we have to do a partial dot product.
# For the moment, do this the explicit, non-optimized way
return (self.conjugate()*x).sum(spaces=spaces)
def s_vdot(self, x):
"""Computes the dot product of 'self' with x.
Parameters
----------
x : Field
x must be defined on the same domain as `self`.
Returns
-------
float or complex
The dot product
"""
if not isinstance(x, Field):
raise TypeError("The dot-partner must be an instance of " +
"the Field class")
if x._domain != self._domain:
raise ValueError("Domain mismatch")
return np.vdot(self._val, x._val)
def norm(self, ord=2):
"""Computes the L2-norm of the field values.
......@@ -357,7 +379,7 @@ class Field(object):
def _contraction_helper(self, op, spaces):
if spaces is None:
return getattr(self._val, op)()
return Field.scalar(getattr(self._val, op)())
spaces = utilities.parse_spaces(spaces, len(self._domain))
......@@ -371,7 +393,7 @@ class Field(object):
# check if the result is scalar or if a result_field must be constr.
if np.isscalar(data):
return data
return Field.scalar(data)
else:
return_domain = tuple(dom
for i, dom in enumerate(self._domain)
......@@ -390,12 +412,21 @@ class Field(object):
Returns
-------
Field or scalar
The result of the summation. If it is carried out over the entire
domain, this is a scalar, otherwise a Field.
Field
The result of the summation.
"""
return self._contraction_helper('sum', spaces)
def s_sum(self):
"""Returns the sum over all entries
Returns
-------
scalar
The result of the summation.
"""
return self._val.sum()
def integrate(self, spaces=None):
"""Integrates over the sub-domains given by `spaces`.
......@@ -410,9 +441,8 @@ class Field(object):
Returns
-------
Field or scalar
The result of the integration. If it is carried out over the
entire domain, this is a scalar, otherwise a Field.
Field
The result of the integration.
"""
swgt = self.scalar_weight(spaces)
if swgt is not None:
......@@ -422,6 +452,23 @@ class Field(object):
tmp = self.weight(1, spaces=spaces)
return tmp.sum(spaces)
def s_integrate(self):
"""Integrates over the Field.
Integration is performed by summing over `self` multiplied by its
volume factors.
Returns
-------
Scalar
The result of the integration.
"""
swgt = self.scalar_weight()
if swgt is not None:
return self.s_sum()*swgt
tmp = self.weight(1)
return tmp.s_sum()
def prod(self, spaces=None):
"""Computes the product over the sub-domains given by `spaces`.
......@@ -434,18 +481,26 @@ class Field(object):
Returns
-------
Field or scalar
The result of the product. If it is carried out over the entire
domain, this is a scalar, otherwise a Field.
Field
The result of the product.
"""
return self._contraction_helper('prod', spaces)
def s_prod(self):
return self._val.prod()
def all(self, spaces=None):
return self._contraction_helper('all', spaces)
def s_all(self):
return self._val.all()
def any(self, spaces=None):
return self._contraction_helper('any', spaces)
def s_any(self):
return self._val.any()
# def min(self, spaces=None):
# """Determines the minimum over the sub-domains given by `spaces`.
#
......@@ -457,9 +512,8 @@ class Field(object):
#
# Returns
# -------
# Field or scalar
# The result of the operation. If it is carried out over the entire
# domain, this is a scalar, otherwise a Field.
# Field
# The result of the operation.
# """
# return self._contraction_helper('min', spaces)
#
......@@ -474,9 +528,8 @@ class Field(object):
#
# Returns
# -------
# Field or scalar
# The result of the operation. If it is carried out over the entire
# domain, this is a scalar, otherwise a Field.
# Field
# The result of the operation.
# """
# return self._contraction_helper('max', spaces)
......@@ -494,9 +547,8 @@ class Field(object):
Returns
-------
Field or scalar
The result of the operation. If it is carried out over the entire
domain, this is a scalar, otherwise a Field.
Field
The result of the operation.
"""
if self.scalar_weight(spaces) is not None:
return self._contraction_helper('mean', spaces)
......@@ -505,6 +557,19 @@ class Field(object):
tmp = self.weight(1, spaces)
return tmp.sum(spaces)*(1./tmp.total_volume(spaces))
def s_mean(self):
"""Determines the field mean
``x.s_mean()`` is equivalent to
``x.s_integrate()/x.total_volume()``.
Returns
-------
scalar
The result of the operation.
"""
return self.s_integrate()/self.total_volume()