Commit f827fba2 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'NIFTy_4' into fix_nonlinear_cg

parents 57a67f20 9b934fba
# custom
# never store the git version file
git_version.py
# custom
setup.cfg
.idea
.DS_Store
......
......@@ -82,9 +82,9 @@ author = u'Theo Steininger / Martin Reinecke'
# built documents.
#
# The short X.Y version.
version = u'3.9'
version = u'4.0'
# The full version, including alpha/beta/rc tags.
release = u'3.9.0'
release = u'4.0.0'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
......
......@@ -16,52 +16,27 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from .. import Field, exp
from ..operators.diagonal_operator import DiagonalOperator
from ..field import Field, exp
from ..minimization.energy import Energy
# TODO Take only residual_sample_list as argument
from ..operators.diagonal_operator import DiagonalOperator
class NoiseEnergy(Energy):
def __init__(self, position, d, xi, D, t, ht, Instrument,
nonlinearity, alpha, q, Distributor, samples=3,
xi_sample_list=None, inverter=None):
super(NoiseEnergy, self).__init__(position=position)
self.xi = xi
self.D = D
self.d = d
self.N = DiagonalOperator(diagonal=exp(self.position))
self.t = t
self.samples = samples
self.ht = ht
self.Instrument = Instrument
self.nonlinearity = nonlinearity
def __init__(self, position, alpha, q, res_sample_list):
super(NoiseEnergy, self).__init__(position)
self.N = DiagonalOperator(diagonal=exp(self.position))
self.alpha = alpha
self.q = q
self.Distributor = Distributor
self.power = self.Distributor(exp(0.5 * self.t))
if xi_sample_list is None:
if samples is None or samples == 0:
xi_sample_list = [xi]
else:
xi_sample_list = [D.draw_sample() + xi
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
A = Distributor(exp(.5*self.t))
alpha_field = Field(self.position.domain, val=alpha)
q_field = Field(self.position.domain, val=q)
self.res_sample_list = res_sample_list
self._gradient = None
for sample in self.xi_sample_list:
map_s = self.ht(A * sample)
residual = self.d - \
self.Instrument(self.nonlinearity(map_s))
lh = .5 * residual.vdot(self.N.inverse_times(residual))
grad = -.5 * self.N.inverse_times(residual.conjugate()*residual)
for s in self.res_sample_list:
lh = .5 * s.vdot(self.N.inverse_times(s))
grad = -.5 * self.N.inverse_times(s.conjugate()*s)
if self._gradient is None:
self._value = lh
self._gradient = grad.copy()
......@@ -69,20 +44,19 @@ class NoiseEnergy(Energy):
self._value += lh
self._gradient += grad
self._value *= 1. / len(self.xi_sample_list)
expmpos = exp(-position)
self._value *= 1./len(self.res_sample_list)
self._value += .5 * self.position.sum()
self._value += (self.alpha - 1.).vdot(self.position) + \
self.q.vdot(exp(-self.position))
self._value += (alpha_field-1.).vdot(self.position) + \
q_field.vdot(expmpos)
self._gradient *= 1. / len(self.xi_sample_list)
self._gradient += (self.alpha-0.5) - self.q*(exp(-self.position))
self._gradient *= 1./len(self.res_sample_list)
self._gradient += (alpha_field-0.5) - q_field*expmpos
self._gradient.lock()
def at(self, position):
return self.__class__(
position, self.d, self.xi, self.D, self.t, self.ht,
self.Instrument, self.nonlinearity, self.alpha, self.q,
self.Distributor, xi_sample_list=self.xi_sample_list,
samples=self.samples, inverter=self.inverter)
return self.__class__(position, self.alpha, self.q,
self.res_sample_list)
@property
def value(self):
......
......@@ -81,11 +81,13 @@ class DescentMinimizer(Minimizer):
# compute a step length that reduces energy.value sufficiently
try:
new_energy = self.line_searcher.perform_line_search(
new_energy, success = self.line_searcher.perform_line_search(
energy=energy, pk=self.get_descent_direction(energy),
f_k_minus_1=f_k_minus_1)
except ValueError:
return energy, controller.ERROR
if not success:
self.reset()
f_k_minus_1 = energy.value
......@@ -103,6 +105,9 @@ class DescentMinimizer(Minimizer):
if status != controller.CONTINUE:
return energy, status
def reset(self):
pass
@abc.abstractmethod
def get_descent_direction(self, energy):
""" Calculates the next descent direction.
......
......@@ -109,3 +109,19 @@ class Energy(NiftyMetaBase()):
curvature of the potential at the given `position`.
"""
raise NotImplementedError
def longest_step(self, dir):
"""Returns the longest allowed step size along `dir`
Parameters
----------
dir : Field
the search direction
Returns
-------
float or None
the longest allowed step when starting from `self.position` along
`dir`. If None, the step size is not limited.
"""
return None
......@@ -54,5 +54,7 @@ class LineSearch(NiftyMetaBase()):
-------
Energy
The new Energy object on the new position.
bool
whether the line search was considered successful or not
"""
raise NotImplementedError
......@@ -41,7 +41,7 @@ class LineSearchStrongWolfe(LineSearch):
Parameter for curvature condition rule. (Default: 0.9)
max_step_size : float
Maximum step allowed in to be made in the descent direction.
(Default: 1000000000)
(Default: 1e30)
max_iterations : int, optional
Maximum number of iterations performed by the line search algorithm.
(Default: 100)
......@@ -51,7 +51,7 @@ class LineSearchStrongWolfe(LineSearch):
"""
def __init__(self, preferred_initial_step_size=None, c1=1e-4, c2=0.9,
max_step_size=1000000000, max_iterations=100,
max_step_size=1e30, max_iterations=100,
max_zoom_iterations=100):
super(LineSearchStrongWolfe, self).__init__(
......@@ -85,19 +85,26 @@ class LineSearchStrongWolfe(LineSearch):
-------
Energy
The new Energy object on the new position.
bool
whether the line search was considered successful or not
"""
le_0 = LineEnergy(0., energy, pk, 0.)
maxstepsize = energy.longest_step(pk)
if maxstepsize is None:
maxstepsize = self.max_step_size
maxstepsize = min(maxstepsize, self.max_step_size)
# initialize the zero phis
old_phi_0 = f_k_minus_1
phi_0 = le_0.value
phiprime_0 = le_0.directional_derivative
if phiprime_0 == 0:
dobj.mprint("Directional derivative is zero; assuming convergence")
return energy
return energy, False
if phiprime_0 > 0:
dobj.mprint("Error: search direction is not a descent direction")
raise ValueError("search direction must be a descent direction")
return energy, False
# set alphas
alpha0 = 0.
......@@ -112,49 +119,44 @@ class LineSearchStrongWolfe(LineSearch):
alpha1 = 1.0
else:
alpha1 = 1.0/pk.norm()
alpha1 = min(alpha1, 0.99*maxstepsize)
# start the minimization loop
iteration_number = 0
while iteration_number < self.max_iterations:
iteration_number += 1
if alpha1 == 0:
result_energy = le_0.energy
break
return le_0.energy, False
le_alpha1 = le_0.at(alpha1)
phi_alpha1 = le_alpha1.value
if (phi_alpha1 > phi_0 + self.c1*alpha1*phiprime_0) or \
((phi_alpha1 >= phi_alpha0) and (iteration_number > 1)):
le_star = self._zoom(alpha0, alpha1, phi_0, phiprime_0,
phi_alpha0, phiprime_alpha0, phi_alpha1,
le_0)
result_energy = le_star.energy
break
return self._zoom(alpha0, alpha1, phi_0, phiprime_0,
phi_alpha0, phiprime_alpha0, phi_alpha1,
le_0)
phiprime_alpha1 = le_alpha1.directional_derivative
if abs(phiprime_alpha1) <= -self.c2*phiprime_0:
result_energy = le_alpha1.energy
break
return le_alpha1.energy, True
if phiprime_alpha1 >= 0:
le_star = self._zoom(alpha1, alpha0, phi_0, phiprime_0,
phi_alpha1, phiprime_alpha1, phi_alpha0,
le_0)
result_energy = le_star.energy
break
return self._zoom(alpha1, alpha0, phi_0, phiprime_0,
phi_alpha1, phiprime_alpha1, phi_alpha0,
le_0)
# update alphas
alpha0, alpha1 = alpha1, min(2*alpha1, self.max_step_size)
if alpha1 == self.max_step_size:
return le_alpha1.energy
alpha0, alpha1 = alpha1, min(2*alpha1, maxstepsize)
if alpha1 == maxstepsize:
dobj.mprint("max step size reached")
return le_alpha1.energy, False
phi_alpha0 = phi_alpha1
phiprime_alpha0 = phiprime_alpha1
else:
dobj.mprint("max iterations reached")
return le_alpha1.energy
return result_energy
dobj.mprint("max iterations reached")
return le_alpha1.energy, False
def _zoom(self, alpha_lo, alpha_hi, phi_0, phiprime_0,
phi_lo, phiprime_lo, phi_hi, le_0):
......@@ -238,7 +240,7 @@ class LineSearchStrongWolfe(LineSearch):
phiprime_alphaj = le_alphaj.directional_derivative
# If the second Wolfe condition is met, return the result
if abs(phiprime_alphaj) <= -self.c2*phiprime_0:
return le_alphaj
return le_alphaj.energy, True
# If not, check the sign of the slope
if phiprime_alphaj*delta_alpha >= 0:
alpha_recent, phi_recent = alpha_hi, phi_hi
......@@ -251,7 +253,7 @@ class LineSearchStrongWolfe(LineSearch):
else:
dobj.mprint("The line search algorithm (zoom) did not converge.")
return le_alphaj
return le_alphaj.energy, False
def _cubicmin(self, a, fa, fpa, b, fb, c, fc):
"""Estimating the minimum with cubic interpolation.
......
......@@ -63,8 +63,10 @@ class NonlinearCG(Minimizer):
while True:
grad_old = energy.gradient
f_k = energy.value
energy = self._line_searcher.perform_line_search(energy, p,
f_k_minus_1)
energy, success = self._line_searcher.perform_line_search(
energy, p, f_k_minus_1)
if not success:
return energy, controller.ERROR
f_k_minus_1 = f_k
status = self._controller.check(energy)
if status != controller.CONTINUE:
......
......@@ -16,7 +16,7 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from __future__ import division, print_function
from __future__ import division
from .minimizer import Minimizer
from ..field import Field
from .. import dobj
......@@ -97,9 +97,9 @@ class ScipyMinimizer(Minimizer):
status = self._controller.check(hlp._energy)
return hlp._energy, self._controller.check(hlp._energy)
if not r.success:
print("Problem in Scipy minimization:", r.message)
dobj.mprint("Problem in Scipy minimization:", r.message)
else:
print("Problem in Scipy minimization")
dobj.mprint("Problem in Scipy minimization")
return hlp._energy, self._controller.ERROR
......
......@@ -49,6 +49,9 @@ class VL_BFGS(DescentMinimizer):
self._information_store = None
return super(VL_BFGS, self).__call__(energy)
def reset(self):
self._information_store = None
def get_descent_direction(self, energy):
x = energy.position
gradient = energy.gradient
......
......@@ -138,11 +138,10 @@ class LinearOperator(NiftyMetaBase()):
other = self._toOperator(other, self.domain)
return SumOperator.make([self, other], [False, True])
# MR FIXME: this might be more complicated ...
def __rsub__(self, other):
from .sum_operator import SumOperator
other = self._toOperator(other, self.domain)
return SumOperator.make(other, self, [False, True])
return SumOperator.make([other, self], [False, True])
@abc.abstractproperty
def capability(self):
......
......@@ -198,6 +198,18 @@ def plot(f, **kwargs):
if not isinstance(label, list):
label = [label]
linewidth = _get_kw("linewidth", None, **kwargs)
if linewidth is None:
linewidth = [1.] * len(f)
if not isinstance(linewidth, list):
linewidth = [linewidth]
alpha = _get_kw("alpha", None, **kwargs)
if alpha is None:
alpha = [None] * len(f)
if not isinstance(alpha, list):
alpha = [alpha]
dom = dom[0]
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
......@@ -216,7 +228,8 @@ def plot(f, **kwargs):
xcoord = np.arange(npoints, dtype=np.float64)*dist
for i, fld in enumerate(f):
ycoord = fld.to_global_data()
plt.plot(xcoord, ycoord, label=label[i])
plt.plot(xcoord, ycoord, label=label[i],
linewidth=linewidth[i], alpha=alpha[i])
_limit_xy(**kwargs)
if label != ([None]*len(f)):
plt.legend()
......@@ -249,7 +262,8 @@ def plot(f, **kwargs):
xcoord = dom.k_lengths
for i, fld in enumerate(f):
ycoord = fld.to_global_data()
plt.plot(xcoord, ycoord, label=label[i])
plt.plot(xcoord, ycoord, label=label[i],
linewidth=linewidth[i], alpha=alpha[i])
_limit_xy(**kwargs)
if label != ([None]*len(f)):
plt.legend()
......
......@@ -2,4 +2,13 @@
# 1) we don't load dependencies by storing it in __init__.py
# 2) we can import it in setup.py for the same reason
# 3) we can import it into your module module
__version__ = '3.9.0'
__version__ = '4.0.0'
def gitversion():
try:
from .git_version import gitversion
except ImportError:
return "unknown"
return gitversion
......@@ -18,6 +18,17 @@
from setuptools import setup, find_packages
def write_version():
import subprocess
p = subprocess.Popen(["git", "describe", "--dirty", "--tags"],
stdout=subprocess.PIPE)
res = p.communicate()[0].strip().decode('utf-8')
with open("nifty4/git_version.py", "w") as file:
file.write('gitversion = "{}"\n'.format(res))
write_version()
exec(open('nifty4/version.py').read())
setup(name="nifty4",
......
......@@ -73,7 +73,7 @@ class Noise_Energy_Tests(unittest.TestCase):
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
D = ift.library.NonlinearWienerFilterEnergy(
C = ift.library.NonlinearWienerFilterEnergy(
position=xi,
d=d,
Instrument=R,
......@@ -84,10 +84,10 @@ class Noise_Energy_Tests(unittest.TestCase):
S=S,
inverter=inverter).curvature
energy0 = ift.library.NoiseEnergy(
position=eta0, d=d, xi=xi, D=D, t=tau, Instrument=R,
alpha=alpha, q=q, Distributor=Dist, nonlinearity=f,
ht=ht, samples=10)
res_sample_list = [d - R(f(ht(C.draw_sample() + xi)))
for _ in range(10)]
energy0 = ift.library.NoiseEnergy(eta0, alpha, q, res_sample_list)
energy1 = energy0.at(eta1)
a = (energy1.value - energy0.value) / eps
......
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