Commit 88872f64 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge nifty2go

parents cb10770a f89a5b17
Pipeline #19152 canceled with stage
......@@ -88,11 +88,11 @@ if __name__ == "__main__":
d_data = d.val.real
ift.plotting.plot(d.real, name="data.pdf")
IC1 = ift.DefaultIterationController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
IC1 = ift.GradientNormController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
minimizer1 = ift.RelaxedNewton(IC1)
IC2 = ift.DefaultIterationController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
IC2 = ift.GradientNormController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
IC3 = ift.DefaultIterationController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
IC3 = ift.GradientNormController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
minimizer3 = ift.SteepestDescent(IC3)
# Set starting position
......@@ -107,14 +107,14 @@ if __name__ == "__main__":
S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
# Initialize non-linear Wiener Filter energy
ICI = ift.DefaultIterationController(verbose=False,iteration_limit=500,tol_abs_gradnorm=0.1)
ICI = ift.GradientNormController(verbose=False,iteration_limit=500,tol_abs_gradnorm=0.1)
map_inverter = ift.ConjugateGradient(controller=ICI)
map_energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0, inverter=map_inverter)
# Solve the Wiener Filter analytically
D0 = map_energy.curvature
m0 = D0.inverse_times(j)
# Initialize power energy with updated parameters
ICI2 = ift.DefaultIterationController(name="powI",verbose=True,iteration_limit=200,tol_abs_gradnorm=1e-5)
ICI2 = ift.GradientNormController(name="powI",verbose=True,iteration_limit=200,tol_abs_gradnorm=1e-5)
power_inverter = ift.ConjugateGradient(controller=ICI2)
power_energy = ift.library.CriticalPowerEnergy(position=t0, m=m0, D=D0,
smoothness_prior=10., samples=3, inverter=power_inverter)
......
......@@ -46,8 +46,8 @@ if __name__ == "__main__":
# Wiener filter
m0 = ift.Field(harmonic_space, val=0.)
ctrl = ift.DefaultIterationController(verbose=False,tol_abs_gradnorm=1)
ctrl2 = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1, name="outer")
ctrl = ift.GradientNormController(verbose=False,tol_abs_gradnorm=1)
ctrl2 = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1, name="outer")
inverter = ift.ConjugateGradient(controller=ctrl)
energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S, inverter=inverter)
minimizer1 = ift.VL_BFGS(controller=ctrl2,max_history_length=20)
......
......@@ -93,7 +93,7 @@ if __name__ == "__main__":
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.DefaultIterationController(verbose=True, tol_custom=1e-3, convergence_level=3)
ctrl = ift.GradientNormController(verbose=True, iteration_limit=100)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic, inverter=inverter)
......
......@@ -45,7 +45,7 @@ if __name__ == "__main__":
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1)
ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic,inverter=inverter)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
......
......@@ -57,7 +57,7 @@ if __name__ == "__main__":
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=1e-2)
ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=1e-2)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic, inverter=inverter)
......
......@@ -75,7 +75,7 @@ if __name__ == "__main__":
# Choosing the minimization strategy
ctrl = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1)
ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
# Setting starting position
m0 = ift.Field(h_space, val=.0)
......
......@@ -19,7 +19,6 @@
from __future__ import division
import abc
from .nifty_meta import NiftyMeta
from future.utils import with_metaclass
......@@ -120,7 +119,7 @@ class DomainObject(with_metaclass(
"There is no generic dim for DomainObject.")
@abc.abstractmethod
def scalar_weight(self):
def scalar_dvol(self):
""" Returns the volume factors of this domain as a floating
point scalar, if the volume factors are all identical, otherwise
returns None.
......@@ -137,10 +136,9 @@ class DomainObject(with_metaclass(
"""
raise NotImplementedError(
"There is no generic scalar_weight method for DomainObject.")
"There is no generic scalar_dvol method for DomainObject.")
@abc.abstractmethod
def weight(self):
def dvol(self):
""" Returns the volume factors of this domain, either as a floating
point scalar (if the volume factors are all identical) or as a
floating point array with a shape of `self.shape`.
......@@ -157,5 +155,4 @@ class DomainObject(with_metaclass(
If called for this abstract class.
"""
raise NotImplementedError(
"There is no generic weight method for DomainObject.")
return self.scalar_dvol()
......@@ -66,7 +66,6 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (object,), {}))):
def __init__(self, position):
super(Energy, self).__init__()
self._cache = {}
self._position = position.copy()
def at(self, position):
......
......@@ -95,10 +95,8 @@ class LineEnergy(object):
"""
return self.__class__(line_position,
self.energy,
self.line_direction,
offset=self.line_position)
return LineEnergy(line_position, self.energy, self.line_direction,
offset=self.line_position)
@property
def value(self):
......
......@@ -18,9 +18,11 @@
def memo(f):
name = id(f)
name = f.__name__
def wrapped_f(self):
if not hasattr(self, "_cache"):
self._cache = {}
try:
return self._cache[name]
except KeyError:
......
......@@ -8,23 +8,21 @@ class QuadraticEnergy(Energy):
position-independent.
"""
def __init__(self, position, A, b, _grad=None, _bnorm=None):
def __init__(self, position, A, b, _grad=None):
super(QuadraticEnergy, self).__init__(position=position)
self._A = A
self._b = b
self._bnorm = _bnorm
if _grad is not None:
self._Ax = _grad + self._b
else:
self._Ax = self._A(self.position)
def at(self, position):
return self.__class__(position=position, A=self._A, b=self._b,
_bnorm=self.norm_b)
return QuadraticEnergy(position=position, A=self._A, b=self._b)
def at_with_grad(self, position, grad):
return self.__class__(position=position, A=self._A, b=self._b,
_grad=grad, _bnorm=self.norm_b)
return QuadraticEnergy(position=position, A=self._A, b=self._b,
_grad=grad)
@property
@memo
......@@ -39,9 +37,3 @@ class QuadraticEnergy(Energy):
@property
def curvature(self):
return self._A
@property
def norm_b(self):
if self._bnorm is None:
self._bnorm = self._b.norm()
return self._bnorm
......@@ -380,15 +380,6 @@ class Field(object):
"""
return self.domain.dim
@property
def total_volume(self):
""" Returns the total volume of all spaces in the domain.
"""
if len(self.domain) == 0:
return 0.
volume_tuple = tuple(sp.total_volume for sp in self.domain)
return reduce(lambda x, y: x * y, volume_tuple)
@property
def real(self):
""" The real part of the field (data is not copied).
......@@ -417,13 +408,13 @@ class Field(object):
def scalar_weight(self, spaces=None):
if np.isscalar(spaces):
return self.domain[spaces].scalar_weight()
return self.domain[spaces].scalar_dvol()
if spaces is None:
spaces = range(len(self.domain))
res = 1.
for i in spaces:
tmp = self.domain[i].scalar_weight()
tmp = self.domain[i].scalar_dvol()
if tmp is None:
return None
res *= tmp
......@@ -459,7 +450,7 @@ class Field(object):
fct = 1.
for ind in spaces:
wgt = self.domain[ind].weight()
wgt = self.domain[ind].dvol()
if np.isscalar(wgt):
fct *= wgt
else:
......
......@@ -20,8 +20,5 @@ from ..domain_object import DomainObject
class FieldType(DomainObject):
def scalar_weight(self):
return 1.
def weight(self):
def scalar_dvol(self):
return 1.
......@@ -26,7 +26,6 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
"""
def __init__(self, R, N, S, d, position, inverter, fft4exp=None, **kwargs):
self._cache = {}
self.R = R
self.N = N
self.S = S
......
......@@ -17,8 +17,7 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from .line_searching import *
from .iteration_controller import IterationController
from .default_iteration_controller import DefaultIterationController
from .iteration_controlling import *
from .minimizer import Minimizer
from .conjugate_gradient import ConjugateGradient
from .nonlinear_cg import NonlinearCG
......
......@@ -73,7 +73,6 @@ class ConjugateGradient(Minimizer):
if status != controller.CONTINUE:
return energy, status
norm_b = energy.norm_b
r = energy.gradient
if preconditioner is not None:
d = preconditioner(r)
......@@ -111,9 +110,7 @@ class ConjugateGradient(Minimizer):
if gamma == 0:
return energy, controller.CONVERGED
status = self._controller.check(energy,
custom_measure=np.sqrt(gamma) /
norm_b)
status = self._controller.check(energy)
if status != controller.CONTINUE:
return energy, status
......
from .iteration_controller import IterationController
from .gradient_norm_controller import GradientNormController
......@@ -20,28 +20,27 @@ from __future__ import print_function
from .iteration_controller import IterationController
class DefaultIterationController(IterationController):
class GradientNormController(IterationController):
def __init__(self, tol_abs_gradnorm=None, tol_rel_gradnorm=None,
tol_custom=None, convergence_level=1, iteration_limit=None,
convergence_level=1, iteration_limit=None,
name=None, verbose=None):
super(DefaultIterationController, self).__init__()
super(GradientNormController, self).__init__()
self._tol_abs_gradnorm = tol_abs_gradnorm
self._tol_rel_gradnorm = tol_rel_gradnorm
self._tol_custom = tol_custom
self._convergence_level = convergence_level
self._iteration_limit = iteration_limit
self._name = name
self._verbose = verbose
def start(self, energy, custom_measure=None):
def start(self, energy):
self._itcount = -1
self._ccount = 0
if self._tol_rel_gradnorm is not None:
self._tol_rel_gradnorm_now = self._tol_rel_gradnorm \
* energy.gradient_norm
return self.check(energy, custom_measure)
return self.check(energy)
def check(self, energy, custom_measure=None):
def check(self, energy):
self._itcount += 1
inclvl = False
......@@ -51,9 +50,6 @@ class DefaultIterationController(IterationController):
if self._tol_rel_gradnorm is not None:
if energy.gradient_norm <= self._tol_rel_gradnorm_now:
inclvl = True
if self._tol_custom is not None and custom_measure is not None:
if custom_measure <= self._tol_custom:
inclvl = True
if inclvl:
self._ccount += 1
else:
......@@ -67,8 +63,6 @@ class DefaultIterationController(IterationController):
msg += " Iteration #" + str(self._itcount)
msg += " energy=" + str(energy.value)
msg += " gradnorm=" + str(energy.gradient_norm)
if custom_measure is not None:
msg += " custom=" + str(custom_measure)
msg += " clvl=" + str(self._ccount)
print(msg)
# self.logger.info(msg)
......
......@@ -18,7 +18,7 @@
from builtins import range
import abc
from ..nifty_meta import NiftyMeta
from ...nifty_meta import NiftyMeta
from future.utils import with_metaclass
......
......@@ -62,9 +62,9 @@ class RGRGTransformation(Transformation):
# BUT: the pixel volumes of the domain and codomain are different.
# Hence, in order to produce the same scalar product, power==1.
if self.codomain.harmonic:
fct = self.domain.weight()
fct = self.domain.scalar_dvol()
else:
fct = 1./(self.codomain.weight()*self.domain.dim)
fct = 1./(self.codomain.scalar_dvol()*self.domain.dim)
# Perform the transformation
if issubclass(val.dtype.type, np.complexfloating):
......
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