Skip to content

GitLab

N
NIFTy
Commit 92ff0e59 authored by Martin Reinecke's avatar Martin Reinecke
Browse files
Options

byebye models

parent 5572575f
Hide whitespace changes
Inline Side-by-side
Showing with 23 additions and 839 deletions
nifty5/__init__.py
View file @ 92ff0e59
......@@ -18,14 +18,6 @@ from .field import Field
from .nonlinearities import Exponential, Linear, PositiveTanh, Tanh
from .models.constant import Constant
from .models.linear_model import LinearModel
from .models.local_nonlinearity import (LocalModel, PointwiseExponential,
PointwisePositiveTanh, PointwiseTanh)
from .models.model import Model
from .models.multi_model import MultiModel
from .models.variable import Variable
from .operators.central_zero_padder import CentralZeroPadder
from .operators.diagonal_operator import DiagonalOperator
from .operators.dof_distributor import DOFDistributor
......@@ -43,14 +35,12 @@ from .operators.inversion_enabler import InversionEnabler
from .operators.laplace_operator import LaplaceOperator
from .operators.linear_operator import LinearOperator
from .operators.mask_operator import MaskOperator
from .operators.multi_adaptor import MultiAdaptor
from .operators.null_operator import NullOperator
from .operators.power_distributor import PowerDistributor
from .operators.qht_operator import QHTOperator
from .operators.sampling_enabler import SamplingEnabler
from .operators.sandwich_operator import SandwichOperator
from .operators.scaling_operator import ScalingOperator
from .operators.selection_operator import SelectionOperator
from .operators.slope_operator import SlopeOperator
from .operators.smoothness_operator import SmoothnessOperator
from .operators.symmetrizing_operator import SymmetrizingOperator
......@@ -81,14 +71,14 @@ from .minimization.energy_sum import EnergySum
from .sugar import *
from .plotting.plot import plot, plot_finish
from .library.amplitude_model import make_amplitude_model, AmplitudeModel
from .library.amplitude_model import AmplitudeModel
from .library.gaussian_energy import GaussianEnergy
from .library.los_response import LOSResponse
from .library.point_sources import PointSources
#from .library.point_sources import PointSources
from .library.poissonian_energy import PoissonianEnergy
from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.correlated_fields import (make_correlated_field,
make_mf_correlated_field)
#from .library.correlated_fields import (make_correlated_field,
# make_mf_correlated_field)
from .library.bernoulli_energy import BernoulliEnergy
from . import extra
......
nifty5/extra/energy_and_model_tests.py
View file @ 92ff0e59
......@@ -20,60 +20,12 @@ from __future__ import absolute_import, division, print_function
from ..compat import *
import numpy as np
from ..sugar import from_random
from ..minimization.energy import Energy
from ..models.model import Model
from ..linearization import Linearization
__all__ = ["check_value_gradient_consistency",
"check_value_gradient_metric_consistency",
"check_value_gradient_metric_consistency2",
"check_value_gradient_consistency2"]
"check_value_gradient_metric_consistency"]
def _get_acceptable_model(M):
val = M.value
if not np.isfinite(val.sum()):
raise ValueError('Initial Model value must be finite')
dir = from_random("normal", M.position.domain)
dirder = M.jacobian(dir)
if dirder.norm() == 0:
dir = dir * val.norm() * 1e-5
else:
dir = dir * val.norm() * (1e-5/dirder.norm())
# Find a step length that leads to a "reasonable" Model
for i in range(50):
try:
M2 = M.at(M.position+dir)
if np.isfinite(M2.value.sum()) and abs(M2.value.sum()) < 1e20:
break
except FloatingPointError:
pass
dir = dir*0.5
else:
raise ValueError("could not find a reasonable initial step")
return M2
def _get_acceptable_energy(E):
val = E.value
if not np.isfinite(val):
raise ValueError('Initial Energy must be finite')
dir = from_random("normal", E.position.domain)
dirder = E.gradient.vdot(dir)
dir = dir * (np.abs(val)/np.abs(dirder)*1e-5)
# Find a step length that leads to a "reasonable" energy
for i in range(50):
try:
E2 = E.at(E.position+dir)
if np.isfinite(E2.value) and abs(E2.value) < 1e20:
break
except FloatingPointError:
pass
dir = dir*0.5
else:
raise ValueError("could not find a reasonable initial step")
return E2
def _get_acceptable_location(op, loc, lin):
val = lin.val
if not np.isfinite(val.sum()):
......@@ -99,39 +51,7 @@ def _get_acceptable_location(op, loc, lin):
return loc2, lin2
def check_value_gradient_consistency(E, tol=1e-8, ntries=100):
for _ in range(ntries):
if isinstance(E, Energy):
E2 = _get_acceptable_energy(E)
else:
E2 = _get_acceptable_model(E)
val = E.value
dir = E2.position - E.position
Enext = E2
dirnorm = dir.norm()
for i in range(50):
Emid = E.at(E.position + 0.5*dir)
if isinstance(E, Energy):
dirder = Emid.gradient.vdot(dir)/dirnorm
else:
dirder = Emid.jacobian(dir)/dirnorm
numgrad = (E2.value-val)/dirnorm
if isinstance(E, Model):
xtol = tol * dirder.norm() / np.sqrt(dirder.size)
if (abs(numgrad-dirder) <= xtol).all():
break
else:
xtol = tol*Emid.gradient_norm
if abs(numgrad-dirder) <= xtol:
break
dir = dir*0.5
dirnorm *= 0.5
E2 = Emid
else:
raise ValueError("gradient and value seem inconsistent")
E = Enext
def check_value_gradient_consistency2(op, loc, tol=1e-8, ntries=100):
def check_value_gradient_consistency(op, loc, tol=1e-8, ntries=100):
for _ in range(ntries):
lin = op(Linearization.make_var(loc))
loc2, lin2 = _get_acceptable_location(op, loc, lin)
......@@ -154,7 +74,9 @@ def check_value_gradient_consistency2(op, loc, tol=1e-8, ntries=100):
else:
raise ValueError("gradient and value seem inconsistent")
loc = locnext
def check_value_gradient_metric_consistency2(op, loc, tol=1e-8, ntries=100):
def check_value_gradient_metric_consistency(op, loc, tol=1e-8, ntries=100):
for _ in range(ntries):
lin = op(Linearization.make_var(loc))
loc2, lin2 = _get_acceptable_location(op, loc, lin)
......@@ -180,27 +102,3 @@ def check_value_gradient_metric_consistency2(op, loc, tol=1e-8, ntries=100):
else:
raise ValueError("gradient and value seem inconsistent")
loc = locnext
def check_value_gradient_metric_consistency(E, tol=1e-8, ntries=100):
if isinstance(E, Model):
raise ValueError('Models have no metric, thus it cannot be tested.')
for _ in range(ntries):
E2 = _get_acceptable_energy(E)
val = E.value
dir = E2.position - E.position
Enext = E2
dirnorm = dir.norm()
for i in range(50):
Emid = E.at(E.position + 0.5*dir)
dirder = Emid.gradient.vdot(dir)/dirnorm
dgrad = Emid.metric(dir)/dirnorm
xtol = tol*Emid.gradient_norm
if abs((E2.value-val)/dirnorm - dirder) < xtol and \
(abs((E2.gradient-E.gradient)/dirnorm-dgrad) < xtol).all():
break
dir = dir*0.5
dirnorm *= 0.5
E2 = Emid
else:
raise ValueError("gradient, value and metric seem inconsistent")
E = Enext
nifty5/library/amplitude_model.py
View file @ 92ff0e59
......@@ -34,74 +34,6 @@ def _ceps_kernel(dof_space, k, a, k0):
return a**2/(1+(k/(k0*dof_space.bindistances[0]))**2)**2
def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi']):
'''
Method for construction of amplitude model
Computes a smooth power spectrum.
Output lives in PowerSpace.
Parameters
----------
Npixdof : #pix in dof_space
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0
sm, sv : slope_mean = expected exponent of power law (e.g. -4),
slope_variance (default=1)
im, iv : y-intercept_mean, y-intercept_variance of power_slope
'''
from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..models.variable import Variable
from ..models.constant import Constant
from ..models.local_nonlinearity import PointwiseExponential
h_space = s_space.get_default_codomain()
p_space = PowerSpace(h_space)
exp_transform = ExpTransform(p_space, Npixdof)
logk_space = exp_transform.domain[0]
qht = QHTOperator(target=logk_space)
dof_space = qht.domain[0]
param_space = UnstructuredDomain(2)
sym = SymmetrizingOperator(logk_space)
phi_mean = np.array([sm, im])
phi_sig = np.array([sv, iv])
slope = SlopeOperator(param_space, logk_space, phi_sig)
norm_phi_mean = Field.from_global_data(param_space, phi_mean/phi_sig)
fields = {keys[0]: Field.from_random('normal', dof_space),
keys[1]: Field.from_random('normal', param_space)}
position = MultiField.from_dict(fields)
dof_space = position[keys[0]].domain[0]
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
ceps = makeOp(sqrt(cepstrum))
smooth_op = sym * qht * ceps
smooth_spec = smooth_op(Variable(position)[keys[0]])
phi = Variable(position)[keys[1]] + Constant(position, norm_phi_mean)
linear_spec = slope(phi)
loglog_spec = smooth_spec + linear_spec
xlog_ampl = PointwiseExponential(0.5*loglog_spec)
internals = {'loglog_spec': loglog_spec,
'qht': qht,
'ceps': ceps,
'norm_phi_mean': norm_phi_mean}
return exp_transform(xlog_ampl), internals
def create_cepstrum_amplitude_field(domain, cepstrum):
"""Creates a ...
Writes the sum of all modes into the zero-mode.
......@@ -147,6 +79,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
return Field.from_global_data(domain, cepstrum_field)
class AmplitudeModel(Operator):
'''
Computes a smooth power spectrum.
......@@ -172,9 +105,6 @@ class AmplitudeModel(Operator):
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
from ..models.variable import Variable
from ..models.constant import Constant
from ..models.local_nonlinearity import PointwiseExponential
h_space = s_space.get_default_codomain()
p_space = PowerSpace(h_space)
......@@ -192,11 +122,7 @@ class AmplitudeModel(Operator):
self._norm_phi_mean = Field.from_global_data(param_space, phi_mean/phi_sig)
self._domain = MultiDomain.make({keys[0]: dof_space, keys[1]: param_space})
# fields = {keys[0]: Field.from_random('normal', dof_space),
# keys[1]: Field.from_random('normal', param_space)}
# position = MultiField.from_dict(fields)
# dof_space = position[keys[0]].domain[0]
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
......@@ -204,13 +130,6 @@ class AmplitudeModel(Operator):
self._smooth_op = sym * qht * ceps
self._keys = tuple(keys)
# smooth_spec = smooth_op(Variable(position)[keys[0]])
# phi = Variable(position)[keys[1]] + Constant(position, norm_phi_mean)
# linear_spec = slope(phi)
# loglog_spec = smooth_spec + linear_spec
# xlog_ampl = PointwiseExponential(0.5*loglog_spec)
def __call__(self, x):
smooth_spec = self._smooth_op(x[self._keys[0]])
phi = x[self._keys[1]] + self._norm_phi_mean
......
nifty5/models/binary_helpers.py deleted 100644 → 0
View file @ 5572575f
# 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 ..multi.multi_field import MultiField
from ..sugar import makeOp
from .model import Model
def _joint_position(model1, model2):
a = model1.position.to_dict()
b = model2.position.to_dict()
# Note: In python >3.5 one could do {**a, **b}
ab = a
ab.update(b)
return MultiField.from_dict(ab)
class ScalarMul(Model):
"""Class representing a model multiplied by a scalar factor."""
def __init__(self, factor, model):
super(ScalarMul, self).__init__(model.position)
# TODO -> floating
if not isinstance(factor, (float, int)):
raise TypeError
self._model = model
self._factor = factor
self._value = self._factor * self._model.value
self._jacobian = self._factor * self._model.jacobian
def at(self, position):
return self.__class__(self._factor, self._model.at(position))
class Add(Model):
"""Class representing the sum of two models."""
def __init__(self, position, model1, model2):
super(Add, self).__init__(position)
self._model1 = model1.at(position)
self._model2 = model2.at(position)
self._value = self._model1.value + self._model2.value
self._jacobian = self._model1.jacobian + self._model2.jacobian
@staticmethod
def make(model1, model2):
"""Build the sum of two models.
Parameters
----------
model1: Model
First model.
model2: Model
Second model
"""
position = _joint_position(model1, model2)
return Add(position, model1, model2)
def at(self, position):
return self.__class__(position, self._model1, self._model2)
class Mul(Model):
"""Class representing the pointwise product of two models."""
def __init__(self, position, model1, model2):
super(Mul, self).__init__(position)
self._model1 = model1.at(position)
self._model2 = model2.at(position)
self._value = self._model1.value * self._model2.value
self._jacobian = (makeOp(self._model1.value) * self._model2.jacobian +
makeOp(self._model2.value) * self._model1.jacobian)
@staticmethod
def make(model1, model2):
"""Build the pointwise product of two models.
Parameters
----------
model1: Model
First model.
model2: Model
Second model
"""
position = _joint_position(model1, model2)
return Mul(position, model1, model2)
def at(self, position):
return self.__class__(position, self._model1, self._model2)
nifty5/models/constant.py deleted 100644 → 0
View file @ 5572575f
# 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 ..operators.null_operator import NullOperator
from .model import Model
class Constant(Model):
"""A model with a constant (multi-)field as value.
Parameters
----------
position : Field or MultiField
The current position in parameter space.
constant : Field
The value of the model.
Notes
-----
Since there is no model-function associated:
- Position has no influence on value.
- The Jacobian is a null matrix.
"""
def __init__(self, position, constant):
super(Constant, self).__init__(position)
self._constant = constant
self._value = self._constant
self._jacobian = NullOperator(position.domain, constant.domain)
def at(self, position):
return self.__class__(position, self._constant)
nifty5/models/linear_model.py deleted 100644 → 0
View file @ 5572575f
# 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 ..operators.selection_operator import SelectionOperator
from .model import Model
class LinearModel(Model):
def __init__(self, inp, lin_op):
"""Computes lin_op(inp) where lin_op is a Linear Operator
Parameters
----------
inp : Model
lin_op : LinearOperator
linear function to be applied to model
Returns
-------
Model with linear Operator applied:
- Model.value = LinOp (inp.value) [key-wise]
- Jacobian = LinOp * inp.jacobian
"""
from ..operators.linear_operator import LinearOperator
super(LinearModel, self).__init__(inp.position)
if not isinstance(lin_op, LinearOperator):
raise TypeError("needs a LinearOperator as input")
self._lin_op = lin_op
self._inp = inp
# MR FIXME: what does this do?
if isinstance(self._lin_op, SelectionOperator):
self._lin_op = SelectionOperator(self._inp.value.domain,
self._lin_op._key)
self._value = self._lin_op(self._inp.value)
self._jacobian = self._lin_op*self._inp.jacobian
def at(self, position):
return self.__class__(self._inp.at(position), self._lin_op)
nifty5/models/local_nonlinearity.py deleted 100644 → 0
View file @ 5572575f
# 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 ..nonlinearities import Exponential, PositiveTanh, Tanh
from ..sugar import makeOp
from .model import Model
class LocalModel(Model):
def __init__(self, inp, nonlinearity):
"""
Computes nonlinearity(inp)
- LocalModel.value = nonlinearity(value) (pointwise)
- LocalModel.jacobian = Outer Product of Jacobians
Parameters
----------
inp : Model
The model for which the nonlinarity will be applied.
nonlinearity: Function
The nonlinearity to be applied to model.
"""
super(LocalModel, self).__init__(inp.position)
self._inp = inp
self._nonlinearity = nonlinearity
self._value = nonlinearity(self._inp.value)
d_inner = self._inp.jacobian
d_outer = makeOp(self._nonlinearity.derivative(self._inp.value))
self._jacobian = d_outer * d_inner
def at(self, position):
return self.__class__(self._inp.at(position), self._nonlinearity)
def PointwiseExponential(inp):
return LocalModel(inp, Exponential())
def PointwiseTanh(inp):
return LocalModel(inp, Tanh())
def PointwisePositiveTanh(inp):
return LocalModel(inp, PositiveTanh())
nifty5/models/model.py deleted 100644 → 0
View file @ 5572575f
# 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 ..field import Field
from ..multi.multi_field import MultiField
from ..operators.selection_operator import SelectionOperator
from ..sugar import makeOp
from ..utilities import NiftyMetaBase
# MR FIXME: it seems that the code expects that the position of a Model is
# always a MultiField ... should we document this?
class Model(NiftyMetaBase()):
"""
The Model object is an implementation of a * which knows:
- position in parameter space. (Field, MultiField)
- value according to its model function A. A(position)
- Jacobian of the model function at the current position.
Parameters
----------
position : Field, MultiField
The input parameter of the model
Notes
-----
An instance of the model class knows its position, value and Jacobian.
One can 'jump' to a new position with the help of the 'at' method, whereby
one automatically gets the value and Jacobian of the model. The 'at' method
creates a new instance of the class.
"""
def __init__(self, position):
self._position = position
def at(self, position):
""" Returns a new Model object, initialized at `position`.
Parameters
----------
position : Field or MultiField
Location in parameter space for the new Model object.
Returns
-------
Model
Model object at new position.
"""
raise NotImplementedError
@property
def position(self):
"""
Field or MultiField: selected location in parameter space.
The location in parameter space where value and Jacobian are
evaluated.
"""
return self._position
@property
def value(self):
<