Commit 9728ac00 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'operator_spectra' into 'NIFTy_5'

Operator spectra

See merge request !333
parents e6188b31 2147bf7c
Pipeline #62097 canceled with stages
in 37 seconds
......@@ -46,15 +46,16 @@ from .operators.outer_product_operator import OuterProduct
from .operators.simple_linear_operators import (
VdotOperator, ConjugationOperator, Realizer,
FieldAdapter, ducktape, GeometryRemover, NullOperator,
MatrixProductOperator, PartialExtractor)
from .operators.value_inserter import ValueInserter
from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
BernoulliEnergy, StandardHamiltonian, AveragedEnergy)
BernoulliEnergy, StandardHamiltonian, AveragedEnergy, QuadraticFormOperator,
from .operators.convolution_operators import FuncConvolutionOperator
from .probing import probe_with_posterior_samples, probe_diagonal, \
StatCalculator, approximation2endo
from .minimization.line_search import LineSearch
from .minimization.iteration_controllers import (
......@@ -97,6 +98,8 @@ from .logger import logger
from .linearization import Linearization
from .operator_spectrum import operator_spectrum
from . import internal_config
_scheme = internal_config.parallelization_scheme()
if _scheme == "Samples":
# 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
# 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 <>.
# Copyright(C) 2013-2019 Max-Planck-Society
import numpy as np
import scipy.sparse.linalg as ssl
from .domain_tuple import DomainTuple
from .domains.unstructured_domain import UnstructuredDomain
from .field import Field
from .multi_domain import MultiDomain
from .multi_field import MultiField
from .operators.linear_operator import LinearOperator
from .operators.sandwich_operator import SandwichOperator
from .sugar import from_global_data, makeDomain
class _DomRemover(LinearOperator):
"""Operator which transforms between a structured MultiDomain
and an unstructured domain.
domain: MultiDomain
the full input domain of the operator.
The operator converts the full domain of its input domain to an
def __init__(self, domain):
self._domain = makeDomain(domain)
if isinstance(self._domain, MultiDomain):
self._size_array = np.array([0] +
[d.size for d in domain.values()])
self._size_array = np.array([0, domain.size])
np.cumsum(self._size_array, out=self._size_array)
target = UnstructuredDomain(self._size_array[-1])
self._target = makeDomain(target)
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
x = x.to_global_data()
if isinstance(self._domain, DomainTuple):
res = x.ravel() if mode == self.TIMES else x.reshape(
res = np.empty( if mode == self.TIMES else {}
for ii, (kk, dd) in enumerate(self.domain.items()):
i0, i1 = self._size_array[ii:ii + 2]
if mode == self.TIMES:
res[i0:i1] = x[kk].ravel()
res[kk] = x[i0:i1].reshape(dd.shape)
return from_global_data(self._tgt(mode), res)
def _check_float_dtype(fld):
if isinstance(fld, MultiField):
dts = [ff.local_data.dtype for ff in fld.values()]
elif isinstance(fld, Field):
dts = [fld.local_data.dtype]
raise TypeError
for dt in dts:
if not np.issubdtype(dt, np.float64):
raise TypeError('Operator supports only floating point dtypes')
def operator_spectrum(A, k, hermitian, which='LM', tol=0):
Find k eigenvalues and eigenvectors of the endomorphism A.
A : LinearOperator
Operator of which eigenvalues shall be computed.
k : int
The number of eigenvalues and eigenvectors desired. `k` must be
smaller than N-1. It is not possible to compute all eigenvectors of a
hermitian: bool
Specifies whether A is hermitian or not.
which : str, ['LM' | 'SM' | 'LR' | 'SR' | 'LI' | 'SI'], optional
Which `k` eigenvectors and eigenvalues to find:
'LM' : largest magnitude
'SM' : smallest magnitude
'LR' : largest real part
'SR' : smallest real part
'LI' : largest imaginary part
'SI' : smallest imaginary part
tol : float, optional
Relative accuracy for eigenvalues (stopping criterion)
The default value of 0 implies machine precision.
w : ndarray
Array of k eigenvalues.
When the requested convergence is not obtained.
The currently converged eigenvalues and eigenvectors can be found
as ``eigenvalues`` and ``eigenvectors`` attributes of the exception
if not isinstance(A, LinearOperator):
raise TypeError('Operator needs to be linear.')
if A.domain is not
raise TypeError('Operator needs to be endomorphism.')
size = A.domain.size
Ar = SandwichOperator.make(_DomRemover(A.domain).adjoint, A)
M = ssl.LinearOperator(
matvec=lambda x: Ar(from_global_data(Ar.domain, x)).to_global_data())
f = ssl.eigsh if hermitian else ssl.eigs
eigs = f(M, k=k, tol=tol, return_eigenvectors=False, which=which)
return np.flip(np.sort(eigs), axis=0)
......@@ -319,9 +319,11 @@ class StandardHamiltonian(EnergyOperator):
def __init__(self, lh, ic_samp=None):
def __init__(self, lh, ic_samp=None, _c_inp=None):
self._lh = lh
self._prior = GaussianEnergy(domain=lh.domain)
if _c_inp is not None:
_, self._prior = self._prior.simplify_for_constant_input(_c_inp)
self._ic_samp = ic_samp
self._domain = lh.domain
......@@ -338,9 +340,13 @@ class StandardHamiltonian(EnergyOperator):
def __repr__(self):
subs = 'Likelihood:\n{}'.format(utilities.indent(self._lh.__repr__()))
subs += '\nPrior: Quadratic{}'.format(self._lh.domain.keys())
subs += '\nPrior:\n{}'.format(self._prior)
return 'StandardHamiltonian:\n' + utilities.indent(subs)
def _simplify_for_constant_input_nontrivial(self, c_inp):
out, lh1 = self._lh.simplify_for_constant_input(c_inp)
return out, StandardHamiltonian(lh1, self._ic_samp, _c_inp=c_inp)
class AveragedEnergy(EnergyOperator):
"""Averages an energy over samples.
......@@ -23,6 +23,7 @@ from ..minimization.iteration_controllers import IterationController
from ..minimization.quadratic_energy import QuadraticEnergy
from ..sugar import full
from .endomorphic_operator import EndomorphicOperator
from .linear_operator import LinearOperator
class InversionEnabler(EndomorphicOperator):
......@@ -47,6 +48,12 @@ class InversionEnabler(EndomorphicOperator):
def __init__(self, op, iteration_controller, approximation=None):
# isinstance(op, EndomorphicOperator) does not suffice since op can be
# a ChainOperator
if not isinstance(op, LinearOperator):
raise TypeError('Operator needs to be linear.')
if op.domain is not
raise TypeError('Operator needs to be endomorphic.')
self._op = op
self._ic = iteration_controller
self._approximation = approximation
......@@ -42,15 +42,20 @@ class SamplingEnabler(EndomorphicOperator):
operator, which supports the operation modes that the operator doesn't
have. It is used as a preconditioner during the iterative inversion,
to accelerate convergence.
start_from_zero : boolean
If true, the conjugate gradient algorithm starts from a field filled
with zeros. Otherwise, it starts from a prior samples. Default is
def __init__(self, likelihood, prior, iteration_controller,
self._op = likelihood + prior
approximation=None, start_from_zero=False):
self._likelihood = likelihood
self._prior = prior
self._ic = iteration_controller
self._approximation = approximation
self._start_from_zero = bool(start_from_zero)
self._op = likelihood + prior
self._domain = self._op.domain
self._capability = self._op.capability
......@@ -60,11 +65,15 @@ class SamplingEnabler(EndomorphicOperator):
except NotImplementedError:
if not from_inverse:
raise ValueError("from_inverse must be True here")
s = self._prior.draw_sample(from_inverse=True)
sp = self._prior(s)
nj = self._likelihood.draw_sample()
energy = QuadraticEnergy(s, self._op, sp + nj,
_grad=self._likelihood(s) - nj)
if self._start_from_zero:
b = self._op.draw_sample()
energy = QuadraticEnergy(0*b, self._op, b)
s = self._prior.draw_sample(from_inverse=True)
sp = self._prior(s)
nj = self._likelihood.draw_sample()
energy = QuadraticEnergy(s, self._op, sp + nj,
_grad=self._likelihood(s) - nj)
inverter = ConjugateGradient(self._ic)
if self._approximation is not None:
energy, convergence = inverter(
......@@ -347,8 +347,12 @@ def _plot2D(f, ax, **kwargs):
if len(dom) == 2:
if (not isinstance(dom[1], RGSpace)) or len(dom[1].shape) != 1:
raise TypeError("need 1D RGSpace as second domain")
rgb = _rgb_data(f.to_global_data())
have_rgb = True
if dom[1].shape[0] == 1:
from .sugar import from_global_data
f = from_global_data(f.domain[0], f.to_global_data()[..., 0])
rgb = _rgb_data(f.to_global_data())
have_rgb = True
foo = kwargs.pop("norm", None)
norm = {} if foo is None else {'norm': foo}
......@@ -477,6 +481,17 @@ class Plot(object):
alpha: float or list of floats
Transparency value.
from .multi_field import MultiField
if isinstance(f, MultiField):
for kk in f.domain.keys():
mykwargs = kwargs.copy()
if 'title' in kwargs:
mykwargs['title'] = "{} {}".format(kk, kwargs['title'])
mykwargs['title'] = "{}".format(kk)
Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment